Bug Summary

File:liboctave/array/CMatrix.cc
Location:line 1884, column 3
Description:Undefined or garbage value returned to caller

Annotated Source Code

1// Matrix manipulations.
2/*
3
4Copyright (C) 1994-2013 John W. Eaton
5Copyright (C) 2008-2009 Jaroslav Hajek
6Copyright (C) 2009 VZLU Prague, a.s.
7
8This file is part of Octave.
9
10Octave is free software; you can redistribute it and/or modify it
11under the terms of the GNU General Public License as published by the
12Free Software Foundation; either version 3 of the License, or (at your
13option) any later version.
14
15Octave is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18for more details.
19
20You should have received a copy of the GNU General Public License
21along with Octave; see the file COPYING. If not, see
22<http://www.gnu.org/licenses/>.
23
24*/
25
26#ifdef HAVE_CONFIG_H1
27#include <config.h>
28#endif
29
30#include <cfloat>
31
32#include <iostream>
33#include <vector>
34
35// FIXME
36#include <sys/types.h>
37
38#include "Array-util.h"
39#include "CMatrix.h"
40#include "CmplxAEPBAL.h"
41#include "CmplxCHOL.h"
42#include "CmplxSCHUR.h"
43#include "CmplxSVD.h"
44#include "DET.h"
45#include "f77-fcn.h"
46#include "functor.h"
47#include "lo-error.h"
48#include "lo-ieee.h"
49#include "lo-mappers.h"
50#include "lo-utils.h"
51#include "mx-base.h"
52#include "mx-cm-dm.h"
53#include "mx-cm-s.h"
54#include "mx-dm-cm.h"
55#include "mx-inlines.cc"
56#include "mx-op-defs.h"
57#include "oct-cmplx.h"
58#include "oct-fftw.h"
59#include "oct-locbuf.h"
60#include "oct-norm.h"
61
62// Fortran functions we call.
63
64extern "C"
65{
66 F77_RET_Tint
67 F77_FUNC (xilaenv, XILAENV)xilaenv_ (const octave_idx_type&,
68 F77_CONST_CHAR_ARG_DECLconst char *,
69 F77_CONST_CHAR_ARG_DECLconst char *,
70 const octave_idx_type&, const octave_idx_type&,
71 const octave_idx_type&, const octave_idx_type&,
72 octave_idx_type&
73 F77_CHAR_ARG_LEN_DECL, long
74 F77_CHAR_ARG_LEN_DECL, long);
75
76 F77_RET_Tint
77 F77_FUNC (zgebal, ZGEBAL)zgebal_ (F77_CONST_CHAR_ARG_DECLconst char *,
78 const octave_idx_type&, Complex*,
79 const octave_idx_type&, octave_idx_type&,
80 octave_idx_type&, double*, octave_idx_type&
81 F77_CHAR_ARG_LEN_DECL, long);
82
83 F77_RET_Tint
84 F77_FUNC (dgebak, DGEBAK)dgebak_ (F77_CONST_CHAR_ARG_DECLconst char *,
85 F77_CONST_CHAR_ARG_DECLconst char *,
86 const octave_idx_type&, const octave_idx_type&,
87 const octave_idx_type&, double*,
88 const octave_idx_type&, double*,
89 const octave_idx_type&, octave_idx_type&
90 F77_CHAR_ARG_LEN_DECL, long
91 F77_CHAR_ARG_LEN_DECL, long);
92
93 F77_RET_Tint
94 F77_FUNC (zgemm, ZGEMM)zgemm_ (F77_CONST_CHAR_ARG_DECLconst char *,
95 F77_CONST_CHAR_ARG_DECLconst char *,
96 const octave_idx_type&, const octave_idx_type&,
97 const octave_idx_type&, const Complex&,
98 const Complex*, const octave_idx_type&,
99 const Complex*, const octave_idx_type&,
100 const Complex&, Complex*, const octave_idx_type&
101 F77_CHAR_ARG_LEN_DECL, long
102 F77_CHAR_ARG_LEN_DECL, long);
103
104 F77_RET_Tint
105 F77_FUNC (zgemv, ZGEMV)zgemv_ (F77_CONST_CHAR_ARG_DECLconst char *,
106 const octave_idx_type&, const octave_idx_type&,
107 const Complex&, const Complex*,
108 const octave_idx_type&, const Complex*,
109 const octave_idx_type&, const Complex&,
110 Complex*, const octave_idx_type&
111 F77_CHAR_ARG_LEN_DECL, long);
112
113 F77_RET_Tint
114 F77_FUNC (xzdotu, XZDOTU)xzdotu_ (const octave_idx_type&, const Complex*,
115 const octave_idx_type&, const Complex*,
116 const octave_idx_type&, Complex&);
117
118 F77_RET_Tint
119 F77_FUNC (xzdotc, XZDOTC)xzdotc_ (const octave_idx_type&, const Complex*,
120 const octave_idx_type&, const Complex*,
121 const octave_idx_type&, Complex&);
122
123 F77_RET_Tint
124 F77_FUNC (zsyrk, ZSYRK)zsyrk_ (F77_CONST_CHAR_ARG_DECLconst char *,
125 F77_CONST_CHAR_ARG_DECLconst char *,
126 const octave_idx_type&, const octave_idx_type&,
127 const Complex&, const Complex*,
128 const octave_idx_type&, const Complex&,
129 Complex*, const octave_idx_type&
130 F77_CHAR_ARG_LEN_DECL, long
131 F77_CHAR_ARG_LEN_DECL, long);
132
133 F77_RET_Tint
134 F77_FUNC (zherk, ZHERK)zherk_ (F77_CONST_CHAR_ARG_DECLconst char *,
135 F77_CONST_CHAR_ARG_DECLconst char *,
136 const octave_idx_type&, const octave_idx_type&,
137 const double&, const Complex*,
138 const octave_idx_type&, const double&, Complex*,
139 const octave_idx_type&
140 F77_CHAR_ARG_LEN_DECL, long
141 F77_CHAR_ARG_LEN_DECL, long);
142
143 F77_RET_Tint
144 F77_FUNC (zgetrf, ZGETRF)zgetrf_ (const octave_idx_type&, const octave_idx_type&,
145 Complex*, const octave_idx_type&,
146 octave_idx_type*, octave_idx_type&);
147
148 F77_RET_Tint
149 F77_FUNC (zgetrs, ZGETRS)zgetrs_ (F77_CONST_CHAR_ARG_DECLconst char *,
150 const octave_idx_type&, const octave_idx_type&,
151 Complex*, const octave_idx_type&,
152 const octave_idx_type*, Complex*,
153 const octave_idx_type&, octave_idx_type&
154 F77_CHAR_ARG_LEN_DECL, long);
155
156 F77_RET_Tint
157 F77_FUNC (zgetri, ZGETRI)zgetri_ (const octave_idx_type&, Complex*,
158 const octave_idx_type&, const octave_idx_type*,
159 Complex*, const octave_idx_type&,
160 octave_idx_type&);
161
162 F77_RET_Tint
163 F77_FUNC (zgecon, ZGECON)zgecon_ (F77_CONST_CHAR_ARG_DECLconst char *,
164 const octave_idx_type&, Complex*,
165 const octave_idx_type&, const double&, double&,
166 Complex*, double*, octave_idx_type&
167 F77_CHAR_ARG_LEN_DECL, long);
168
169 F77_RET_Tint
170 F77_FUNC (zgelsy, ZGELSY)zgelsy_ (const octave_idx_type&, const octave_idx_type&,
171 const octave_idx_type&, Complex*,
172 const octave_idx_type&, Complex*,
173 const octave_idx_type&, octave_idx_type*,
174 double&, octave_idx_type&, Complex*,
175 const octave_idx_type&, double*,
176 octave_idx_type&);
177
178 F77_RET_Tint
179 F77_FUNC (zgelsd, ZGELSD)zgelsd_ (const octave_idx_type&, const octave_idx_type&,
180 const octave_idx_type&, Complex*,
181 const octave_idx_type&, Complex*,
182 const octave_idx_type&, double*, double&,
183 octave_idx_type&, Complex*,
184 const octave_idx_type&, double*,
185 octave_idx_type*, octave_idx_type&);
186
187 F77_RET_Tint
188 F77_FUNC (zpotrf, ZPOTRF)zpotrf_ (F77_CONST_CHAR_ARG_DECLconst char *,
189 const octave_idx_type&, Complex*,
190 const octave_idx_type&, octave_idx_type&
191 F77_CHAR_ARG_LEN_DECL, long);
192
193 F77_RET_Tint
194 F77_FUNC (zpocon, ZPOCON)zpocon_ (F77_CONST_CHAR_ARG_DECLconst char *,
195 const octave_idx_type&, Complex*,
196 const octave_idx_type&, const double&,
197 double&, Complex*, double*, octave_idx_type&
198 F77_CHAR_ARG_LEN_DECL, long);
199
200 F77_RET_Tint
201 F77_FUNC (zpotrs, ZPOTRS)zpotrs_ (F77_CONST_CHAR_ARG_DECLconst char *,
202 const octave_idx_type&, const octave_idx_type&,
203 const Complex*, const octave_idx_type&, Complex*,
204 const octave_idx_type&, octave_idx_type&
205 F77_CHAR_ARG_LEN_DECL, long);
206
207 F77_RET_Tint
208 F77_FUNC (ztrtri, ZTRTRI)ztrtri_ (F77_CONST_CHAR_ARG_DECLconst char *,
209 F77_CONST_CHAR_ARG_DECLconst char *,
210 const octave_idx_type&, const Complex*,
211 const octave_idx_type&, octave_idx_type&
212 F77_CHAR_ARG_LEN_DECL, long
213 F77_CHAR_ARG_LEN_DECL, long);
214
215 F77_RET_Tint
216 F77_FUNC (ztrcon, ZTRCON)ztrcon_ (F77_CONST_CHAR_ARG_DECLconst char *,
217 F77_CONST_CHAR_ARG_DECLconst char *,
218 F77_CONST_CHAR_ARG_DECLconst char *,
219 const octave_idx_type&, const Complex*,
220 const octave_idx_type&, double&,
221 Complex*, double*, octave_idx_type&
222 F77_CHAR_ARG_LEN_DECL, long
223 F77_CHAR_ARG_LEN_DECL, long
224 F77_CHAR_ARG_LEN_DECL, long);
225
226 F77_RET_Tint
227 F77_FUNC (ztrtrs, ZTRTRS)ztrtrs_ (F77_CONST_CHAR_ARG_DECLconst char *,
228 F77_CONST_CHAR_ARG_DECLconst char *,
229 F77_CONST_CHAR_ARG_DECLconst char *,
230 const octave_idx_type&, const octave_idx_type&,
231 const Complex*, const octave_idx_type&, Complex*,
232 const octave_idx_type&, octave_idx_type&
233 F77_CHAR_ARG_LEN_DECL, long
234 F77_CHAR_ARG_LEN_DECL, long
235 F77_CHAR_ARG_LEN_DECL, long);
236
237 F77_RET_Tint
238 F77_FUNC (zlartg, ZLARTG)zlartg_ (const Complex&, const Complex&, double&,
239 Complex&, Complex&);
240
241 F77_RET_Tint
242 F77_FUNC (ztrsyl, ZTRSYL)ztrsyl_ (F77_CONST_CHAR_ARG_DECLconst char *,
243 F77_CONST_CHAR_ARG_DECLconst char *,
244 const octave_idx_type&, const octave_idx_type&,
245 const octave_idx_type&, const Complex*,
246 const octave_idx_type&, const Complex*,
247 const octave_idx_type&, const Complex*,
248 const octave_idx_type&, double&, octave_idx_type&
249 F77_CHAR_ARG_LEN_DECL, long
250 F77_CHAR_ARG_LEN_DECL, long);
251
252 F77_RET_Tint
253 F77_FUNC (xzlange, XZLANGE)xzlange_ (F77_CONST_CHAR_ARG_DECLconst char *,
254 const octave_idx_type&, const octave_idx_type&,
255 const Complex*, const octave_idx_type&,
256 double*, double&
257 F77_CHAR_ARG_LEN_DECL, long);
258}
259
260static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
261
262// Complex Matrix class
263
264ComplexMatrix::ComplexMatrix (const Matrix& a)
265 : MArray<Complex> (a)
266{
267}
268
269ComplexMatrix::ComplexMatrix (const RowVector& rv)
270 : MArray<Complex> (rv)
271{
272}
273
274ComplexMatrix::ComplexMatrix (const ColumnVector& cv)
275 : MArray<Complex> (cv)
276{
277}
278
279ComplexMatrix::ComplexMatrix (const DiagMatrix& a)
280 : MArray<Complex> (a.dims (), 0.0)
281{
282 for (octave_idx_type i = 0; i < a.length (); i++)
283 elem (i, i) = a.elem (i, i);
284}
285
286ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv)
287 : MArray<Complex> (rv)
288{
289}
290
291ComplexMatrix::ComplexMatrix (const ComplexColumnVector& cv)
292 : MArray<Complex> (cv)
293{
294}
295
296ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a)
297 : MArray<Complex> (a.dims (), 0.0)
298{
299 for (octave_idx_type i = 0; i < a.length (); i++)
300 elem (i, i) = a.elem (i, i);
301}
302
303// FIXME: could we use a templated mixed-type copy function here?
304
305ComplexMatrix::ComplexMatrix (const boolMatrix& a)
306 : MArray<Complex> (a)
307{
308}
309
310ComplexMatrix::ComplexMatrix (const charMatrix& a)
311 : MArray<Complex> (a.dims (), 0.0)
312{
313 for (octave_idx_type i = 0; i < a.rows (); i++)
314 for (octave_idx_type j = 0; j < a.cols (); j++)
315 elem (i, j) = static_cast<unsigned char> (a.elem (i, j));
316}
317
318ComplexMatrix::ComplexMatrix (const Matrix& re, const Matrix& im)
319 : MArray<Complex> (re.dims ())
320{
321 if (im.rows () != rows () || im.cols () != cols ())
322 (*current_liboctave_error_handler) ("complex: internal error");
323
324 octave_idx_type nel = numel ();
325 for (octave_idx_type i = 0; i < nel; i++)
326 xelem (i) = Complex (re(i), im(i));
327}
328
329bool
330ComplexMatrix::operator == (const ComplexMatrix& a) const
331{
332 if (rows () != a.rows () || cols () != a.cols ())
333 return false;
334
335 return mx_inline_equal (length (), data (), a.data ());
336}
337
338bool
339ComplexMatrix::operator != (const ComplexMatrix& a) const
340{
341 return !(*this == a);
342}
343
344bool
345ComplexMatrix::is_hermitian (void) const
346{
347 octave_idx_type nr = rows ();
348 octave_idx_type nc = cols ();
349
350 if (is_square () && nr > 0)
351 {
352 for (octave_idx_type i = 0; i < nr; i++)
353 for (octave_idx_type j = i; j < nc; j++)
354 if (elem (i, j) != conj (elem (j, i)))
355 return false;
356
357 return true;
358 }
359
360 return false;
361}
362
363// destructive insert/delete/reorder operations
364
365ComplexMatrix&
366ComplexMatrix::insert (const Matrix& a, octave_idx_type r, octave_idx_type c)
367{
368 octave_idx_type a_nr = a.rows ();
369 octave_idx_type a_nc = a.cols ();
370
371 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
372 {
373 (*current_liboctave_error_handler) ("range error for insert");
374 return *this;
375 }
376
377 if (a_nr >0 && a_nc > 0)
378 {
379 make_unique ();
380
381 for (octave_idx_type j = 0; j < a_nc; j++)
382 for (octave_idx_type i = 0; i < a_nr; i++)
383 xelem (r+i, c+j) = a.elem (i, j);
384 }
385
386 return *this;
387}
388
389ComplexMatrix&
390ComplexMatrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c)
391{
392 octave_idx_type a_len = a.length ();
393
394 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
395 {
396 (*current_liboctave_error_handler) ("range error for insert");
397 return *this;
398 }
399
400 if (a_len > 0)
401 {
402 make_unique ();
403
404 for (octave_idx_type i = 0; i < a_len; i++)
405 xelem (r, c+i) = a.elem (i);
406 }
407
408 return *this;
409}
410
411ComplexMatrix&
412ComplexMatrix::insert (const ColumnVector& a,
413 octave_idx_type r, octave_idx_type c)
414{
415 octave_idx_type a_len = a.length ();
416
417 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
418 {
419 (*current_liboctave_error_handler) ("range error for insert");
420 return *this;
421 }
422
423 if (a_len > 0)
424 {
425 make_unique ();
426
427 for (octave_idx_type i = 0; i < a_len; i++)
428 xelem (r+i, c) = a.elem (i);
429 }
430
431 return *this;
432}
433
434ComplexMatrix&
435ComplexMatrix::insert (const DiagMatrix& a,
436 octave_idx_type r, octave_idx_type c)
437{
438 octave_idx_type a_nr = a.rows ();
439 octave_idx_type a_nc = a.cols ();
440
441 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
442 {
443 (*current_liboctave_error_handler) ("range error for insert");
444 return *this;
445 }
446
447 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
448
449 octave_idx_type a_len = a.length ();
450
451 if (a_len > 0)
452 {
453 make_unique ();
454
455 for (octave_idx_type i = 0; i < a_len; i++)
456 xelem (r+i, c+i) = a.elem (i, i);
457 }
458
459 return *this;
460}
461
462ComplexMatrix&
463ComplexMatrix::insert (const ComplexMatrix& a,
464 octave_idx_type r, octave_idx_type c)
465{
466 Array<Complex>::insert (a, r, c);
467 return *this;
468}
469
470ComplexMatrix&
471ComplexMatrix::insert (const ComplexRowVector& a,
472 octave_idx_type r, octave_idx_type c)
473{
474 octave_idx_type a_len = a.length ();
475 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
476 {
477 (*current_liboctave_error_handler) ("range error for insert");
478 return *this;
479 }
480
481 for (octave_idx_type i = 0; i < a_len; i++)
482 elem (r, c+i) = a.elem (i);
483
484 return *this;
485}
486
487ComplexMatrix&
488ComplexMatrix::insert (const ComplexColumnVector& a,
489 octave_idx_type r, octave_idx_type c)
490{
491 octave_idx_type a_len = a.length ();
492
493 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
494 {
495 (*current_liboctave_error_handler) ("range error for insert");
496 return *this;
497 }
498
499 if (a_len > 0)
500 {
501 make_unique ();
502
503 for (octave_idx_type i = 0; i < a_len; i++)
504 xelem (r+i, c) = a.elem (i);
505 }
506
507 return *this;
508}
509
510ComplexMatrix&
511ComplexMatrix::insert (const ComplexDiagMatrix& a,
512 octave_idx_type r, octave_idx_type c)
513{
514 octave_idx_type a_nr = a.rows ();
515 octave_idx_type a_nc = a.cols ();
516
517 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
518 {
519 (*current_liboctave_error_handler) ("range error for insert");
520 return *this;
521 }
522
523 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
524
525 octave_idx_type a_len = a.length ();
526
527 if (a_len > 0)
528 {
529 make_unique ();
530
531 for (octave_idx_type i = 0; i < a_len; i++)
532 xelem (r+i, c+i) = a.elem (i, i);
533 }
534
535 return *this;
536}
537
538ComplexMatrix&
539ComplexMatrix::fill (double val)
540{
541 octave_idx_type nr = rows ();
542 octave_idx_type nc = cols ();
543
544 if (nr > 0 && nc > 0)
545 {
546 make_unique ();
547
548 for (octave_idx_type j = 0; j < nc; j++)
549 for (octave_idx_type i = 0; i < nr; i++)
550 xelem (i, j) = val;
551 }
552
553 return *this;
554}
555
556ComplexMatrix&
557ComplexMatrix::fill (const Complex& val)
558{
559 octave_idx_type nr = rows ();
560 octave_idx_type nc = cols ();
561
562 if (nr > 0 && nc > 0)
563 {
564 make_unique ();
565
566 for (octave_idx_type j = 0; j < nc; j++)
567 for (octave_idx_type i = 0; i < nr; i++)
568 xelem (i, j) = val;
569 }
570
571 return *this;
572}
573
574ComplexMatrix&
575ComplexMatrix::fill (double val, octave_idx_type r1, octave_idx_type c1,
576 octave_idx_type r2, octave_idx_type c2)
577{
578 octave_idx_type nr = rows ();
579 octave_idx_type nc = cols ();
580
581 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
582 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
583 {
584 (*current_liboctave_error_handler) ("range error for fill");
585 return *this;
586 }
587
588 if (r1 > r2) { std::swap (r1, r2); }
589 if (c1 > c2) { std::swap (c1, c2); }
590
591 if (r2 >= r1 && c2 >= c1)
592 {
593 make_unique ();
594
595 for (octave_idx_type j = c1; j <= c2; j++)
596 for (octave_idx_type i = r1; i <= r2; i++)
597 xelem (i, j) = val;
598 }
599
600 return *this;
601}
602
603ComplexMatrix&
604ComplexMatrix::fill (const Complex& val, octave_idx_type r1, octave_idx_type c1,
605 octave_idx_type r2, octave_idx_type c2)
606{
607 octave_idx_type nr = rows ();
608 octave_idx_type nc = cols ();
609
610 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
611 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
612 {
613 (*current_liboctave_error_handler) ("range error for fill");
614 return *this;
615 }
616
617 if (r1 > r2) { std::swap (r1, r2); }
618 if (c1 > c2) { std::swap (c1, c2); }
619
620 if (r2 >= r1 && c2 >=c1)
621 {
622 make_unique ();
623
624 for (octave_idx_type j = c1; j <= c2; j++)
625 for (octave_idx_type i = r1; i <= r2; i++)
626 xelem (i, j) = val;
627 }
628
629 return *this;
630}
631
632ComplexMatrix
633ComplexMatrix::append (const Matrix& a) const
634{
635 octave_idx_type nr = rows ();
636 octave_idx_type nc = cols ();
637 if (nr != a.rows ())
638 {
639 (*current_liboctave_error_handler) ("row dimension mismatch for append");
640 return *this;
641 }
642
643 octave_idx_type nc_insert = nc;
644 ComplexMatrix retval (nr, nc + a.cols ());
645 retval.insert (*this, 0, 0);
646 retval.insert (a, 0, nc_insert);
647 return retval;
648}
649
650ComplexMatrix
651ComplexMatrix::append (const RowVector& a) const
652{
653 octave_idx_type nr = rows ();
654 octave_idx_type nc = cols ();
655 if (nr != 1)
656 {
657 (*current_liboctave_error_handler) ("row dimension mismatch for append");
658 return *this;
659 }
660
661 octave_idx_type nc_insert = nc;
662 ComplexMatrix retval (nr, nc + a.length ());
663 retval.insert (*this, 0, 0);
664 retval.insert (a, 0, nc_insert);
665 return retval;
666}
667
668ComplexMatrix
669ComplexMatrix::append (const ColumnVector& a) const
670{
671 octave_idx_type nr = rows ();
672 octave_idx_type nc = cols ();
673 if (nr != a.length ())
674 {
675 (*current_liboctave_error_handler) ("row dimension mismatch for append");
676 return *this;
677 }
678
679 octave_idx_type nc_insert = nc;
680 ComplexMatrix retval (nr, nc + 1);
681 retval.insert (*this, 0, 0);
682 retval.insert (a, 0, nc_insert);
683 return retval;
684}
685
686ComplexMatrix
687ComplexMatrix::append (const DiagMatrix& a) const
688{
689 octave_idx_type nr = rows ();
690 octave_idx_type nc = cols ();
691 if (nr != a.rows ())
692 {
693 (*current_liboctave_error_handler) ("row dimension mismatch for append");
694 return *this;
695 }
696
697 octave_idx_type nc_insert = nc;
698 ComplexMatrix retval (nr, nc + a.cols ());
699 retval.insert (*this, 0, 0);
700 retval.insert (a, 0, nc_insert);
701 return retval;
702}
703
704ComplexMatrix
705ComplexMatrix::append (const ComplexMatrix& a) const
706{
707 octave_idx_type nr = rows ();
708 octave_idx_type nc = cols ();
709 if (nr != a.rows ())
710 {
711 (*current_liboctave_error_handler) ("row dimension mismatch for append");
712 return *this;
713 }
714
715 octave_idx_type nc_insert = nc;
716 ComplexMatrix retval (nr, nc + a.cols ());
717 retval.insert (*this, 0, 0);
718 retval.insert (a, 0, nc_insert);
719 return retval;
720}
721
722ComplexMatrix
723ComplexMatrix::append (const ComplexRowVector& a) const
724{
725 octave_idx_type nr = rows ();
726 octave_idx_type nc = cols ();
727 if (nr != 1)
728 {
729 (*current_liboctave_error_handler) ("row dimension mismatch for append");
730 return *this;
731 }
732
733 octave_idx_type nc_insert = nc;
734 ComplexMatrix retval (nr, nc + a.length ());
735 retval.insert (*this, 0, 0);
736 retval.insert (a, 0, nc_insert);
737 return retval;
738}
739
740ComplexMatrix
741ComplexMatrix::append (const ComplexColumnVector& a) const
742{
743 octave_idx_type nr = rows ();
744 octave_idx_type nc = cols ();
745 if (nr != a.length ())
746 {
747 (*current_liboctave_error_handler) ("row dimension mismatch for append");
748 return *this;
749 }
750
751 octave_idx_type nc_insert = nc;
752 ComplexMatrix retval (nr, nc + 1);
753 retval.insert (*this, 0, 0);
754 retval.insert (a, 0, nc_insert);
755 return retval;
756}
757
758ComplexMatrix
759ComplexMatrix::append (const ComplexDiagMatrix& a) const
760{
761 octave_idx_type nr = rows ();
762 octave_idx_type nc = cols ();
763 if (nr != a.rows ())
764 {
765 (*current_liboctave_error_handler) ("row dimension mismatch for append");
766 return *this;
767 }
768
769 octave_idx_type nc_insert = nc;
770 ComplexMatrix retval (nr, nc + a.cols ());
771 retval.insert (*this, 0, 0);
772 retval.insert (a, 0, nc_insert);
773 return retval;
774}
775
776ComplexMatrix
777ComplexMatrix::stack (const Matrix& a) const
778{
779 octave_idx_type nr = rows ();
780 octave_idx_type nc = cols ();
781 if (nc != a.cols ())
782 {
783 (*current_liboctave_error_handler)
784 ("column dimension mismatch for stack");
785 return *this;
786 }
787
788 octave_idx_type nr_insert = nr;
789 ComplexMatrix retval (nr + a.rows (), nc);
790 retval.insert (*this, 0, 0);
791 retval.insert (a, nr_insert, 0);
792 return retval;
793}
794
795ComplexMatrix
796ComplexMatrix::stack (const RowVector& a) const
797{
798 octave_idx_type nr = rows ();
799 octave_idx_type nc = cols ();
800 if (nc != a.length ())
801 {
802 (*current_liboctave_error_handler)
803 ("column dimension mismatch for stack");
804 return *this;
805 }
806
807 octave_idx_type nr_insert = nr;
808 ComplexMatrix retval (nr + 1, nc);
809 retval.insert (*this, 0, 0);
810 retval.insert (a, nr_insert, 0);
811 return retval;
812}
813
814ComplexMatrix
815ComplexMatrix::stack (const ColumnVector& a) const
816{
817 octave_idx_type nr = rows ();
818 octave_idx_type nc = cols ();
819 if (nc != 1)
820 {
821 (*current_liboctave_error_handler)
822 ("column dimension mismatch for stack");
823 return *this;
824 }
825
826 octave_idx_type nr_insert = nr;
827 ComplexMatrix retval (nr + a.length (), nc);
828 retval.insert (*this, 0, 0);
829 retval.insert (a, nr_insert, 0);
830 return retval;
831}
832
833ComplexMatrix
834ComplexMatrix::stack (const DiagMatrix& a) const
835{
836 octave_idx_type nr = rows ();
837 octave_idx_type nc = cols ();
838 if (nc != a.cols ())
839 {
840 (*current_liboctave_error_handler)
841 ("column dimension mismatch for stack");
842 return *this;
843 }
844
845 octave_idx_type nr_insert = nr;
846 ComplexMatrix retval (nr + a.rows (), nc);
847 retval.insert (*this, 0, 0);
848 retval.insert (a, nr_insert, 0);
849 return retval;
850}
851
852ComplexMatrix
853ComplexMatrix::stack (const ComplexMatrix& a) const
854{
855 octave_idx_type nr = rows ();
856 octave_idx_type nc = cols ();
857 if (nc != a.cols ())
858 {
859 (*current_liboctave_error_handler)
860 ("column dimension mismatch for stack");
861 return *this;
862 }
863
864 octave_idx_type nr_insert = nr;
865 ComplexMatrix retval (nr + a.rows (), nc);
866 retval.insert (*this, 0, 0);
867 retval.insert (a, nr_insert, 0);
868 return retval;
869}
870
871ComplexMatrix
872ComplexMatrix::stack (const ComplexRowVector& a) const
873{
874 octave_idx_type nr = rows ();
875 octave_idx_type nc = cols ();
876 if (nc != a.length ())
877 {
878 (*current_liboctave_error_handler)
879 ("column dimension mismatch for stack");
880 return *this;
881 }
882
883 octave_idx_type nr_insert = nr;
884 ComplexMatrix retval (nr + 1, nc);
885 retval.insert (*this, 0, 0);
886 retval.insert (a, nr_insert, 0);
887 return retval;
888}
889
890ComplexMatrix
891ComplexMatrix::stack (const ComplexColumnVector& a) const
892{
893 octave_idx_type nr = rows ();
894 octave_idx_type nc = cols ();
895 if (nc != 1)
896 {
897 (*current_liboctave_error_handler)
898 ("column dimension mismatch for stack");
899 return *this;
900 }
901
902 octave_idx_type nr_insert = nr;
903 ComplexMatrix retval (nr + a.length (), nc);
904 retval.insert (*this, 0, 0);
905 retval.insert (a, nr_insert, 0);
906 return retval;
907}
908
909ComplexMatrix
910ComplexMatrix::stack (const ComplexDiagMatrix& a) const
911{
912 octave_idx_type nr = rows ();
913 octave_idx_type nc = cols ();
914 if (nc != a.cols ())
915 {
916 (*current_liboctave_error_handler)
917 ("column dimension mismatch for stack");
918 return *this;
919 }
920
921 octave_idx_type nr_insert = nr;
922 ComplexMatrix retval (nr + a.rows (), nc);
923 retval.insert (*this, 0, 0);
924 retval.insert (a, nr_insert, 0);
925 return retval;
926}
927
928ComplexMatrix
929conj (const ComplexMatrix& a)
930{
931 return do_mx_unary_map<Complex, Complex, std::conj<double> > (a);
932}
933
934// resize is the destructive equivalent for this one
935
936ComplexMatrix
937ComplexMatrix::extract (octave_idx_type r1, octave_idx_type c1,
938 octave_idx_type r2, octave_idx_type c2) const
939{
940 if (r1 > r2) { std::swap (r1, r2); }
941 if (c1 > c2) { std::swap (c1, c2); }
942
943 return index (idx_vector (r1, r2+1), idx_vector (c1, c2+1));
944}
945
946ComplexMatrix
947ComplexMatrix::extract_n (octave_idx_type r1, octave_idx_type c1,
948 octave_idx_type nr, octave_idx_type nc) const
949{
950 return index (idx_vector (r1, r1 + nr), idx_vector (c1, c1 + nc));
951}
952
953// extract row or column i.
954
955ComplexRowVector
956ComplexMatrix::row (octave_idx_type i) const
957{
958 return index (idx_vector (i), idx_vector::colon);
959}
960
961ComplexColumnVector
962ComplexMatrix::column (octave_idx_type i) const
963{
964 return index (idx_vector::colon, idx_vector (i));
965}
966
967ComplexMatrix
968ComplexMatrix::inverse (void) const
969{
970 octave_idx_type info;
971 double rcon;
972 MatrixType mattype (*this);
973 return inverse (mattype, info, rcon, 0, 0);
974}
975
976ComplexMatrix
977ComplexMatrix::inverse (octave_idx_type& info) const
978{
979 double rcon;
980 MatrixType mattype (*this);
981 return inverse (mattype, info, rcon, 0, 0);
982}
983
984ComplexMatrix
985ComplexMatrix::inverse (octave_idx_type& info, double& rcon, int force,
986 int calc_cond) const
987{
988 MatrixType mattype (*this);
989 return inverse (mattype, info, rcon, force, calc_cond);
990}
991
992ComplexMatrix
993ComplexMatrix::inverse (MatrixType &mattype) const
994{
995 octave_idx_type info;
996 double rcon;
997 return inverse (mattype, info, rcon, 0, 0);
998}
999
1000ComplexMatrix
1001ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const
1002{
1003 double rcon;
1004 return inverse (mattype, info, rcon, 0, 0);
1005}
1006
1007ComplexMatrix
1008ComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info,
1009 double& rcon, int force, int calc_cond) const
1010{
1011 ComplexMatrix retval;
1012
1013 octave_idx_type nr = rows ();
1014 octave_idx_type nc = cols ();
1015
1016 if (nr != nc || nr == 0 || nc == 0)
1017 (*current_liboctave_error_handler) ("inverse requires square matrix");
1018 else
1019 {
1020 int typ = mattype.type ();
1021 char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
1022 char udiag = 'N';
1023 retval = *this;
1024 Complex *tmp_data = retval.fortran_vec ();
1025
1026 F77_XFCN (ztrtri, ZTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1027 F77_CONST_CHAR_ARG2 (&udiag, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1028 nr, tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1029 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1030 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1031
1032 // Throw-away extra info LAPACK gives so as to not change output.
1033 rcon = 0.0;
1034 if (info != 0)
1035 info = -1;
1036 else if (calc_cond)
1037 {
1038 octave_idx_type ztrcon_info = 0;
1039 char job = '1';
1040
1041 OCTAVE_LOCAL_BUFFER (Complex, cwork, 2*nr)octave_local_buffer<Complex> _buffer_cwork (2*nr); Complex
*cwork = _buffer_cwork
;
1042 OCTAVE_LOCAL_BUFFER (double, rwork, nr)octave_local_buffer<double> _buffer_rwork (nr); double *
rwork = _buffer_rwork
;
1043
1044 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1045 F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1046 F77_CONST_CHAR_ARG2 (&udiag, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1047 nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1048 cwork, rwork, ztrcon_infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1049 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1050 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1051 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr
, rcon, cwork, rwork, ztrcon_info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1052
1053 if (ztrcon_info != 0)
1054 info = -1;
1055 }
1056
1057 if (info == -1 && ! force)
1058 retval = *this; // Restore matrix contents.
1059 }
1060
1061 return retval;
1062}
1063
1064ComplexMatrix
1065ComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info,
1066 double& rcon, int force, int calc_cond) const
1067{
1068 ComplexMatrix retval;
1069
1070 octave_idx_type nr = rows ();
1071 octave_idx_type nc = cols ();
1072
1073 if (nr != nc)
1074 (*current_liboctave_error_handler) ("inverse requires square matrix");
1075 else
1076 {
1077 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1078 octave_idx_type *pipvt = ipvt.fortran_vec ();
1079
1080 retval = *this;
1081 Complex *tmp_data = retval.fortran_vec ();
1082
1083 Array<Complex> z (dim_vector (1, 1));
1084 octave_idx_type lwork = -1;
1085
1086 // Query the optimum work array size.
1087
1088 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetri_ (nc, tmp_data, nr, pipvt, z.fortran_vec (), lwork
, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1089 z.fortran_vec (), lwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetri_ (nc, tmp_data, nr, pipvt, z.fortran_vec (), lwork
, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1090
1091 lwork = static_cast<octave_idx_type> (std::real (z(0)));
1092 lwork = (lwork < 2 *nc ? 2*nc : lwork);
1093 z.resize (dim_vector (lwork, 1));
1094 Complex *pz = z.fortran_vec ();
1095
1096 info = 0;
1097
1098 // Calculate the norm of the matrix, for later use.
1099 double anorm;
1100 if (calc_cond)
1101 anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0))
1102 .max ();
1103
1104 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrf_ (nc, nc, tmp_data, nr, pipvt, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1105
1106 // Throw-away extra info LAPACK gives so as to not change output.
1107 rcon = 0.0;
1108 if (info != 0)
1109 info = -1;
1110 else if (calc_cond)
1111 {
1112 // Now calculate the condition number for non-singular matrix.
1113 octave_idx_type zgecon_info = 0;
1114 char job = '1';
1115 Array<double> rz (dim_vector (2 * nc, 1));
1116 double *prz = rz.fortran_vec ();
1117 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, zgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1118 nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, zgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1119 rcon, pz, prz, zgecon_infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, zgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1120 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, zgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1121
1122 if (zgecon_info != 0)
1123 info = -1;
1124 }
1125
1126 if (info == -1 && ! force)
1127 retval = *this; // Restore contents.
1128 else
1129 {
1130 octave_idx_type zgetri_info = 0;
1131
1132 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetri_ (nc, tmp_data, nr, pipvt, pz, lwork, zgetri_info)
; octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1133 pz, lwork, zgetri_info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetri_ (nc, tmp_data, nr, pipvt, pz, lwork, zgetri_info)
; octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1134
1135 if (zgetri_info != 0)
1136 info = -1;
1137 }
1138
1139 if (info != 0)
1140 mattype.mark_as_rectangular ();
1141 }
1142
1143 return retval;
1144}
1145
1146ComplexMatrix
1147ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info,
1148 double& rcon, int force, int calc_cond) const
1149{
1150 int typ = mattype.type (false);
1151 ComplexMatrix ret;
1152
1153 if (typ == MatrixType::Unknown)
1154 typ = mattype.type (*this);
1155
1156 if (typ == MatrixType::Upper || typ == MatrixType::Lower)
1157 ret = tinverse (mattype, info, rcon, force, calc_cond);
1158 else
1159 {
1160 if (mattype.is_hermitian ())
1161 {
1162 ComplexCHOL chol (*this, info, calc_cond);
1163 if (info == 0)
1164 {
1165 if (calc_cond)
1166 rcon = chol.rcond ();
1167 else
1168 rcon = 1.0;
1169 ret = chol.inverse ();
1170 }
1171 else
1172 mattype.mark_as_unsymmetric ();
1173 }
1174
1175 if (!mattype.is_hermitian ())
1176 ret = finverse (mattype, info, rcon, force, calc_cond);
1177
1178 if ((mattype.is_hermitian () || calc_cond) && rcon == 0.)
1179 ret = ComplexMatrix (rows (), columns (), Complex (octave_Inf, 0.));
1180 }
1181
1182 return ret;
1183}
1184
1185ComplexMatrix
1186ComplexMatrix::pseudo_inverse (double tol) const
1187{
1188 ComplexMatrix retval;
1189
1190 ComplexSVD result (*this, SVD::economy);
1191
1192 DiagMatrix S = result.singular_values ();
1193 ComplexMatrix U = result.left_singular_matrix ();
1194 ComplexMatrix V = result.right_singular_matrix ();
1195
1196 ColumnVector sigma = S.extract_diag ();
1197
1198 octave_idx_type r = sigma.length () - 1;
1199 octave_idx_type nr = rows ();
1200 octave_idx_type nc = cols ();
1201
1202 if (tol <= 0.0)
1203 {
1204 if (nr > nc)
1205 tol = nr * sigma.elem (0) * std::numeric_limits<double>::epsilon ();
1206 else
1207 tol = nc * sigma.elem (0) * std::numeric_limits<double>::epsilon ();
1208 }
1209
1210 while (r >= 0 && sigma.elem (r) < tol)
1211 r--;
1212
1213 if (r < 0)
1214 retval = ComplexMatrix (nc, nr, 0.0);
1215 else
1216 {
1217 ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
1218 DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
1219 ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
1220 retval = Vr * D * Ur.hermitian ();
1221 }
1222
1223 return retval;
1224}
1225
1226#if defined (HAVE_FFTW1)
1227
1228ComplexMatrix
1229ComplexMatrix::fourier (void) const
1230{
1231 size_t nr = rows ();
1232 size_t nc = cols ();
1233
1234 ComplexMatrix retval (nr, nc);
1235
1236 size_t npts, nsamples;
1237
1238 if (nr == 1 || nc == 1)
1239 {
1240 npts = nr > nc ? nr : nc;
1241 nsamples = 1;
1242 }
1243 else
1244 {
1245 npts = nr;
1246 nsamples = nc;
1247 }
1248
1249 const Complex *in (data ());
1250 Complex *out (retval.fortran_vec ());
1251
1252 octave_fftw::fft (in, out, npts, nsamples);
1253
1254 return retval;
1255}
1256
1257ComplexMatrix
1258ComplexMatrix::ifourier (void) const
1259{
1260 size_t nr = rows ();
1261 size_t nc = cols ();
1262
1263 ComplexMatrix retval (nr, nc);
1264
1265 size_t npts, nsamples;
1266
1267 if (nr == 1 || nc == 1)
1268 {
1269 npts = nr > nc ? nr : nc;
1270 nsamples = 1;
1271 }
1272 else
1273 {
1274 npts = nr;
1275 nsamples = nc;
1276 }
1277
1278 const Complex *in (data ());
1279 Complex *out (retval.fortran_vec ());
1280
1281 octave_fftw::ifft (in, out, npts, nsamples);
1282
1283 return retval;
1284}
1285
1286ComplexMatrix
1287ComplexMatrix::fourier2d (void) const
1288{
1289 dim_vector dv(rows (), cols ());
1290
1291 ComplexMatrix retval (rows (), cols ());
1292 const Complex *in (data ());
1293 Complex *out (retval.fortran_vec ());
1294
1295 octave_fftw::fftNd (in, out, 2, dv);
1296
1297 return retval;
1298}
1299
1300ComplexMatrix
1301ComplexMatrix::ifourier2d (void) const
1302{
1303 dim_vector dv(rows (), cols ());
1304
1305 ComplexMatrix retval (rows (), cols ());
1306 const Complex *in (data ());
1307 Complex *out (retval.fortran_vec ());
1308
1309 octave_fftw::ifftNd (in, out, 2, dv);
1310
1311 return retval;
1312}
1313
1314#else
1315
1316extern "C"
1317{
1318 // Note that the original complex fft routines were not written for
1319 // double complex arguments. They have been modified by adding an
1320 // implicit double precision (a-h,o-z) statement at the beginning of
1321 // each subroutine.
1322
1323 F77_RET_Tint
1324 F77_FUNC (zffti, ZFFTI)zffti_ (const octave_idx_type&, Complex*);
1325
1326 F77_RET_Tint
1327 F77_FUNC (zfftf, ZFFTF)zfftf_ (const octave_idx_type&, Complex*, Complex*);
1328
1329 F77_RET_Tint
1330 F77_FUNC (zfftb, ZFFTB)zfftb_ (const octave_idx_type&, Complex*, Complex*);
1331}
1332
1333ComplexMatrix
1334ComplexMatrix::fourier (void) const
1335{
1336 ComplexMatrix retval;
1337
1338 octave_idx_type nr = rows ();
1339 octave_idx_type nc = cols ();
1340
1341 octave_idx_type npts, nsamples;
1342
1343 if (nr == 1 || nc == 1)
1344 {
1345 npts = nr > nc ? nr : nc;
1346 nsamples = 1;
1347 }
1348 else
1349 {
1350 npts = nr;
1351 nsamples = nc;
1352 }
1353
1354 octave_idx_type nn = 4*npts+15;
1355
1356 Array<Complex> wsave (nn, 1);
1357 Complex *pwsave = wsave.fortran_vec ();
1358
1359 retval = *this;
1360 Complex *tmp_data = retval.fortran_vec ();
1361
1362 F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave);
1363
1364 for (octave_idx_type j = 0; j < nsamples; j++)
1365 {
1366 octave_quit ();
1367
1368 F77_FUNC (zfftf, ZFFTF)zfftf_ (npts, &tmp_data[npts*j], pwsave);
1369 }
1370
1371 return retval;
1372}
1373
1374ComplexMatrix
1375ComplexMatrix::ifourier (void) const
1376{
1377 ComplexMatrix retval;
1378
1379 octave_idx_type nr = rows ();
1380 octave_idx_type nc = cols ();
1381
1382 octave_idx_type npts, nsamples;
1383
1384 if (nr == 1 || nc == 1)
1385 {
1386 npts = nr > nc ? nr : nc;
1387 nsamples = 1;
1388 }
1389 else
1390 {
1391 npts = nr;
1392 nsamples = nc;
1393 }
1394
1395 octave_idx_type nn = 4*npts+15;
1396
1397 Array<Complex> wsave (nn, 1);
1398 Complex *pwsave = wsave.fortran_vec ();
1399
1400 retval = *this;
1401 Complex *tmp_data = retval.fortran_vec ();
1402
1403 F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave);
1404
1405 for (octave_idx_type j = 0; j < nsamples; j++)
1406 {
1407 octave_quit ();
1408
1409 F77_FUNC (zfftb, ZFFTB)zfftb_ (npts, &tmp_data[npts*j], pwsave);
1410 }
1411
1412 for (octave_idx_type j = 0; j < npts*nsamples; j++)
1413 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
1414
1415 return retval;
1416}
1417
1418ComplexMatrix
1419ComplexMatrix::fourier2d (void) const
1420{
1421 ComplexMatrix retval;
1422
1423 octave_idx_type nr = rows ();
1424 octave_idx_type nc = cols ();
1425
1426 octave_idx_type npts, nsamples;
1427
1428 if (nr == 1 || nc == 1)
1429 {
1430 npts = nr > nc ? nr : nc;
1431 nsamples = 1;
1432 }
1433 else
1434 {
1435 npts = nr;
1436 nsamples = nc;
1437 }
1438
1439 octave_idx_type nn = 4*npts+15;
1440
1441 Array<Complex> wsave (nn, 1);
1442 Complex *pwsave = wsave.fortran_vec ();
1443
1444 retval = *this;
1445 Complex *tmp_data = retval.fortran_vec ();
1446
1447 F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave);
1448
1449 for (octave_idx_type j = 0; j < nsamples; j++)
1450 {
1451 octave_quit ();
1452
1453 F77_FUNC (zfftf, ZFFTF)zfftf_ (npts, &tmp_data[npts*j], pwsave);
1454 }
1455
1456 npts = nc;
1457 nsamples = nr;
1458 nn = 4*npts+15;
1459
1460 wsave.resize (dim_vector (nn, 1));
1461 pwsave = wsave.fortran_vec ();
1462
1463 Array<Complex> tmp (npts, 1);
1464 Complex *prow = tmp.fortran_vec ();
1465
1466 F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave);
1467
1468 for (octave_idx_type j = 0; j < nsamples; j++)
1469 {
1470 octave_quit ();
1471
1472 for (octave_idx_type i = 0; i < npts; i++)
1473 prow[i] = tmp_data[i*nr + j];
1474
1475 F77_FUNC (zfftf, ZFFTF)zfftf_ (npts, prow, pwsave);
1476
1477 for (octave_idx_type i = 0; i < npts; i++)
1478 tmp_data[i*nr + j] = prow[i];
1479 }
1480
1481 return retval;
1482}
1483
1484ComplexMatrix
1485ComplexMatrix::ifourier2d (void) const
1486{
1487 ComplexMatrix retval;
1488
1489 octave_idx_type nr = rows ();
1490 octave_idx_type nc = cols ();
1491
1492 octave_idx_type npts, nsamples;
1493
1494 if (nr == 1 || nc == 1)
1495 {
1496 npts = nr > nc ? nr : nc;
1497 nsamples = 1;
1498 }
1499 else
1500 {
1501 npts = nr;
1502 nsamples = nc;
1503 }
1504
1505 octave_idx_type nn = 4*npts+15;
1506
1507 Array<Complex> wsave (nn, 1);
1508 Complex *pwsave = wsave.fortran_vec ();
1509
1510 retval = *this;
1511 Complex *tmp_data = retval.fortran_vec ();
1512
1513 F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave);
1514
1515 for (octave_idx_type j = 0; j < nsamples; j++)
1516 {
1517 octave_quit ();
1518
1519 F77_FUNC (zfftb, ZFFTB)zfftb_ (npts, &tmp_data[npts*j], pwsave);
1520 }
1521
1522 for (octave_idx_type j = 0; j < npts*nsamples; j++)
1523 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
1524
1525 npts = nc;
1526 nsamples = nr;
1527 nn = 4*npts+15;
1528
1529 wsave.resize (dim_vector (nn, 1));
1530 pwsave = wsave.fortran_vec ();
1531
1532 Array<Complex> tmp (npts, 1);
1533 Complex *prow = tmp.fortran_vec ();
1534
1535 F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave);
1536
1537 for (octave_idx_type j = 0; j < nsamples; j++)
1538 {
1539 octave_quit ();
1540
1541 for (octave_idx_type i = 0; i < npts; i++)
1542 prow[i] = tmp_data[i*nr + j];
1543
1544 F77_FUNC (zfftb, ZFFTB)zfftb_ (npts, prow, pwsave);
1545
1546 for (octave_idx_type i = 0; i < npts; i++)
1547 tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts);
1548 }
1549
1550 return retval;
1551}
1552
1553#endif
1554
1555ComplexDET
1556ComplexMatrix::determinant (void) const
1557{
1558 octave_idx_type info;
1559 double rcon;
1560 return determinant (info, rcon, 0);
1561}
1562
1563ComplexDET
1564ComplexMatrix::determinant (octave_idx_type& info) const
1565{
1566 double rcon;
1567 return determinant (info, rcon, 0);
1568}
1569
1570ComplexDET
1571ComplexMatrix::determinant (octave_idx_type& info, double& rcon,
1572 int calc_cond) const
1573{
1574 MatrixType mattype (*this);
1575 return determinant (mattype, info, rcon, calc_cond);
1576}
1577
1578ComplexDET
1579ComplexMatrix::determinant (MatrixType& mattype,
1580 octave_idx_type& info, double& rcon,
1581 int calc_cond) const
1582{
1583 ComplexDET retval (1.0);
1584
1585 info = 0;
1586 rcon = 0.0;
1587
1588 octave_idx_type nr = rows ();
1589 octave_idx_type nc = cols ();
1590
1591 if (nr != nc)
1592 (*current_liboctave_error_handler) ("matrix must be square");
1593 else
1594 {
1595 volatile int typ = mattype.type ();
1596
1597 // Even though the matrix is marked as singular (Rectangular), we may
1598 // still get a useful number from the LU factorization, because it always
1599 // completes.
1600
1601 if (typ == MatrixType::Unknown)
1602 typ = mattype.type (*this);
1603 else if (typ == MatrixType::Rectangular)
1604 typ = MatrixType::Full;
1605
1606 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
1607 {
1608 for (octave_idx_type i = 0; i < nc; i++)
1609 retval *= elem (i,i);
1610 }
1611 else if (typ == MatrixType::Hermitian)
1612 {
1613 ComplexMatrix atmp = *this;
1614 Complex *tmp_data = atmp.fortran_vec ();
1615
1616 double anorm = 0;
1617 if (calc_cond) anorm = xnorm (*this, 1);
1618
1619
1620 char job = 'L';
1621 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1622 tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1623 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1624
1625 if (info != 0)
1626 {
1627 rcon = 0.0;
1628 mattype.mark_as_unsymmetric ();
1629 typ = MatrixType::Full;
1630 }
1631 else
1632 {
1633 Array<Complex> z (dim_vector (2 * nc, 1));
1634 Complex *pz = z.fortran_vec ();
1635 Array<double> rz (dim_vector (nc, 1));
1636 double *prz = rz.fortran_vec ();
1637
1638 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1639 nr, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1640 rcon, pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1641 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1642
1643 if (info != 0)
1644 rcon = 0.0;
1645
1646 for (octave_idx_type i = 0; i < nc; i++)
1647 retval *= atmp (i,i);
1648
1649 retval = retval.square ();
1650 }
1651 }
1652 else if (typ != MatrixType::Full)
1653 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1654
1655 if (typ == MatrixType::Full)
1656 {
1657 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1658 octave_idx_type *pipvt = ipvt.fortran_vec ();
1659
1660 ComplexMatrix atmp = *this;
1661 Complex *tmp_data = atmp.fortran_vec ();
1662
1663 info = 0;
1664
1665 // Calculate the norm of the matrix, for later use.
1666 double anorm = 0;
1667 if (calc_cond) anorm = xnorm (*this, 1);
1668
1669 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrf_ (nr, nr, tmp_data, nr, pipvt, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1670
1671 // Throw-away extra info LAPACK gives so as to not change output.
1672 rcon = 0.0;
1673 if (info != 0)
1674 {
1675 info = -1;
1676 retval = ComplexDET ();
1677 }
1678 else
1679 {
1680 if (calc_cond)
1681 {
1682 // Now calc the condition number for non-singular matrix.
1683 char job = '1';
1684 Array<Complex> z (dim_vector (2 * nc, 1));
1685 Complex *pz = z.fortran_vec ();
1686 Array<double> rz (dim_vector (2 * nc, 1));
1687 double *prz = rz.fortran_vec ();
1688
1689 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1690 nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1691 rcon, pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1692 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1693 }
1694
1695 if (info != 0)
1696 {
1697 info = -1;
1698 retval = ComplexDET ();
1699 }
1700 else
1701 {
1702 for (octave_idx_type i = 0; i < nc; i++)
1703 {
1704 Complex c = atmp(i,i);
1705 retval *= (ipvt(i) != (i+1)) ? -c : c;
1706 }
1707 }
1708 }
1709 }
1710 }
1711
1712 return retval;
1713}
1714
1715double
1716ComplexMatrix::rcond (void) const
1717{
1718 MatrixType mattype (*this);
1719 return rcond (mattype);
1720}
1721
1722double
1723ComplexMatrix::rcond (MatrixType &mattype) const
1724{
1725 double rcon;
1
'rcon' declared without an initial value
1726 octave_idx_type nr = rows ();
1727 octave_idx_type nc = cols ();
1728
1729 if (nr != nc)
2
Taking false branch
1730 (*current_liboctave_error_handler) ("matrix must be square");
1731 else if (nr == 0 || nc == 0)
3
Assuming 'nr' is not equal to 0
4
Assuming 'nc' is not equal to 0
5
Taking false branch
1732 rcon = octave_Inf;
1733 else
1734 {
1735 volatile int typ = mattype.type ();
1736
1737 if (typ == MatrixType::Unknown)
6
Assuming 'typ' is not equal to Unknown
7
Taking false branch
1738 typ = mattype.type (*this);
1739
1740 // Only calculate the condition number for LU/Cholesky
1741 if (typ == MatrixType::Upper)
8
Assuming 'typ' is not equal to Upper
9
Taking false branch
1742 {
1743 const Complex *tmp_data = fortran_vec ();
1744 octave_idx_type info = 0;
1745 char norm = '1';
1746 char uplo = 'U';
1747 char dia = 'N';
1748
1749 Array<Complex> z (dim_vector (2 * nc, 1));
1750 Complex *pz = z.fortran_vec ();
1751 Array<double> rz (dim_vector (nc, 1));
1752 double *prz = rz.fortran_vec ();
1753
1754 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1755 F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1756 F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1757 nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1758 pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1759 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1760 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1761 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1762
1763 if (info != 0)
1764 rcon = 0;
1765 }
1766 else if (typ == MatrixType::Permuted_Upper)
10
Assuming 'typ' is not equal to Permuted_Upper
11
Taking false branch
1767 (*current_liboctave_error_handler)
1768 ("permuted triangular matrix not implemented");
1769 else if (typ == MatrixType::Lower)
12
Assuming 'typ' is not equal to Lower
13
Taking false branch
1770 {
1771 const Complex *tmp_data = fortran_vec ();
1772 octave_idx_type info = 0;
1773 char norm = '1';
1774 char uplo = 'L';
1775 char dia = 'N';
1776
1777 Array<Complex> z (dim_vector (2 * nc, 1));
1778 Complex *pz = z.fortran_vec ();
1779 Array<double> rz (dim_vector (nc, 1));
1780 double *prz = rz.fortran_vec ();
1781
1782 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1783 F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1784 F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1785 nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1786 pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1787 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1788 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1789 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1790
1791 if (info != 0)
1792 rcon = 0.0;
1793 }
1794 else if (typ == MatrixType::Permuted_Lower)
14
Assuming 'typ' is not equal to Permuted_Lower
15
Taking false branch
1795 (*current_liboctave_error_handler)
1796 ("permuted triangular matrix not implemented");
1797 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
16
Assuming 'typ' is not equal to Full
17
Assuming 'typ' is equal to Hermitian
18
Taking true branch
1798 {
1799 double anorm = -1.0;
1800
1801 if (typ == MatrixType::Hermitian)
19
Taking true branch
1802 {
1803 octave_idx_type info = 0;
1804 char job = 'L';
1805
1806 ComplexMatrix atmp = *this;
1807 Complex *tmp_data = atmp.fortran_vec ();
1808
1809 anorm = atmp.abs().sum().
1810 row(static_cast<octave_idx_type>(0)).max();
1811
1812 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1813 tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1814 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1815
1816 if (info != 0)
20
Assuming 'info' is equal to 0
21
Taking false branch
1817 {
1818 rcon = 0.0;
1819
1820 mattype.mark_as_unsymmetric ();
1821 typ = MatrixType::Full;
1822 }
1823 else
1824 {
1825 Array<Complex> z (dim_vector (2 * nc, 1));
1826 Complex *pz = z.fortran_vec ();
1827 Array<double> rz (dim_vector (nc, 1));
1828 double *prz = rz.fortran_vec ();
1829
1830 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1831 nr, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1832 rcon, pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1833 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1834
1835 if (info != 0)
22
Taking false branch
1836 rcon = 0.0;
1837 }
1838 }
1839
1840
1841 if (typ == MatrixType::Full)
23
Taking false branch
1842 {
1843 octave_idx_type info = 0;
1844
1845 ComplexMatrix atmp = *this;
1846 Complex *tmp_data = atmp.fortran_vec ();
1847
1848 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1849 octave_idx_type *pipvt = ipvt.fortran_vec ();
1850
1851 if (anorm < 0.)
1852 anorm = atmp.abs ().sum ().
1853 row(static_cast<octave_idx_type>(0)).max ();
1854
1855 Array<Complex> z (dim_vector (2 * nc, 1));
1856 Complex *pz = z.fortran_vec ();
1857 Array<double> rz (dim_vector (2 * nc, 1));
1858 double *prz = rz.fortran_vec ();
1859
1860 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrf_ (nr, nr, tmp_data, nr, pipvt, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1861
1862 if (info != 0)
1863 {
1864 rcon = 0.0;
1865 mattype.mark_as_rectangular ();
1866 }
1867 else
1868 {
1869 char job = '1';
1870 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1871 nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1872 rcon, pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
1873 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
1874
1875 if (info != 0)
1876 rcon = 0.0;
1877 }
1878 }
1879 }
1880 else
1881 rcon = 0.0;
1882 }
1883
1884 return rcon;
24
Undefined or garbage value returned to caller
1885}
1886
1887ComplexMatrix
1888ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
1889 octave_idx_type& info, double& rcon,
1890 solve_singularity_handler sing_handler,
1891 bool calc_cond, blas_trans_type transt) const
1892{
1893 ComplexMatrix retval;
1894
1895 octave_idx_type nr = rows ();
1896 octave_idx_type nc = cols ();
1897
1898 if (nr != b.rows ())
1899 (*current_liboctave_error_handler)
1900 ("matrix dimension mismatch solution of linear equations");
1901 else if (nr == 0 || nc == 0 || b.cols () == 0)
1902 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1903 else
1904 {
1905 volatile int typ = mattype.type ();
1906
1907 if (typ == MatrixType::Permuted_Upper ||
1908 typ == MatrixType::Upper)
1909 {
1910 octave_idx_type b_nc = b.cols ();
1911 rcon = 1.;
1912 info = 0;
1913
1914 if (typ == MatrixType::Permuted_Upper)
1915 {
1916 (*current_liboctave_error_handler)
1917 ("permuted triangular matrix not implemented");
1918 }
1919 else
1920 {
1921 const Complex *tmp_data = fortran_vec ();
1922
1923 if (calc_cond)
1924 {
1925 char norm = '1';
1926 char uplo = 'U';
1927 char dia = 'N';
1928
1929 Array<Complex> z (dim_vector (2 * nc, 1));
1930 Complex *pz = z.fortran_vec ();
1931 Array<double> rz (dim_vector (nc, 1));
1932 double *prz = rz.fortran_vec ();
1933
1934 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1935 F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1936 F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1937 nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1938 pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1939 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1940 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1941 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1942
1943 if (info != 0)
1944 info = -2;
1945
1946 volatile double rcond_plus_one = rcon + 1.0;
1947
1948 if (rcond_plus_one == 1.0 || xisnan (rcon))
1949 {
1950 info = -2;
1951
1952 if (sing_handler)
1953 sing_handler (rcon);
1954 else
1955 (*current_liboctave_error_handler)
1956 ("matrix singular to machine precision, rcond = %g",
1957 rcon);
1958 }
1959 }
1960
1961 if (info == 0)
1962 {
1963 retval = b;
1964 Complex *result = retval.fortran_vec ();
1965
1966 char uplo = 'U';
1967 char trans = get_blas_char (transt);
1968 char dia = 'N';
1969
1970 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1971 F77_CONST_CHAR_ARG2 (&trans, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1972 F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1973 nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1974 result, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1975 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1976 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
1977 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
1978 }
1979 }
1980 }
1981 else
1982 (*current_liboctave_error_handler) ("incorrect matrix type");
1983 }
1984
1985 return retval;
1986}
1987
1988ComplexMatrix
1989ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
1990 octave_idx_type& info, double& rcon,
1991 solve_singularity_handler sing_handler,
1992 bool calc_cond, blas_trans_type transt) const
1993{
1994 ComplexMatrix retval;
1995
1996 octave_idx_type nr = rows ();
1997 octave_idx_type nc = cols ();
1998
1999 if (nr != b.rows ())
2000 (*current_liboctave_error_handler)
2001 ("matrix dimension mismatch solution of linear equations");
2002 else if (nr == 0 || nc == 0 || b.cols () == 0)
2003 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2004 else
2005 {
2006 volatile int typ = mattype.type ();
2007
2008 if (typ == MatrixType::Permuted_Lower ||
2009 typ == MatrixType::Lower)
2010 {
2011 octave_idx_type b_nc = b.cols ();
2012 rcon = 1.;
2013 info = 0;
2014
2015 if (typ == MatrixType::Permuted_Lower)
2016 {
2017 (*current_liboctave_error_handler)
2018 ("permuted triangular matrix not implemented");
2019 }
2020 else
2021 {
2022 const Complex *tmp_data = fortran_vec ();
2023
2024 if (calc_cond)
2025 {
2026 char norm = '1';
2027 char uplo = 'L';
2028 char dia = 'N';
2029
2030 Array<Complex> z (dim_vector (2 * nc, 1));
2031 Complex *pz = z.fortran_vec ();
2032 Array<double> rz (dim_vector (nc, 1));
2033 double *prz = rz.fortran_vec ();
2034
2035 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2036 F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2037 F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2038 nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2039 pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2040 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2041 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2042 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr
, rcon, pz, prz, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2043
2044 if (info != 0)
2045 info = -2;
2046
2047 volatile double rcond_plus_one = rcon + 1.0;
2048
2049 if (rcond_plus_one == 1.0 || xisnan (rcon))
2050 {
2051 info = -2;
2052
2053 if (sing_handler)
2054 sing_handler (rcon);
2055 else
2056 (*current_liboctave_error_handler)
2057 ("matrix singular to machine precision, rcond = %g",
2058 rcon);
2059 }
2060 }
2061
2062 if (info == 0)
2063 {
2064 retval = b;
2065 Complex *result = retval.fortran_vec ();
2066
2067 char uplo = 'L';
2068 char trans = get_blas_char (transt);
2069 char dia = 'N';
2070
2071 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2072 F77_CONST_CHAR_ARG2 (&trans, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2073 F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2074 nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2075 result, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2076 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2077 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2078 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data
, nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2079 }
2080 }
2081 }
2082 else
2083 (*current_liboctave_error_handler) ("incorrect matrix type");
2084 }
2085
2086 return retval;
2087}
2088
2089ComplexMatrix
2090ComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b,
2091 octave_idx_type& info, double& rcon,
2092 solve_singularity_handler sing_handler,
2093 bool calc_cond) const
2094{
2095 ComplexMatrix retval;
2096
2097 octave_idx_type nr = rows ();
2098 octave_idx_type nc = cols ();
2099
2100
2101 if (nr != nc || nr != b.rows ())
2102 (*current_liboctave_error_handler)
2103 ("matrix dimension mismatch solution of linear equations");
2104 else if (nr == 0 || b.cols () == 0)
2105 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2106 else
2107 {
2108 volatile int typ = mattype.type ();
2109
2110 // Calculate the norm of the matrix, for later use.
2111 double anorm = -1.;
2112
2113 if (typ == MatrixType::Hermitian)
2114 {
2115 info = 0;
2116 char job = 'L';
2117
2118 ComplexMatrix atmp = *this;
2119 Complex *tmp_data = atmp.fortran_vec ();
2120
2121 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
2122
2123 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2124 tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2125 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2126
2127 // Throw-away extra info LAPACK gives so as to not change output.
2128 rcon = 0.0;
2129 if (info != 0)
2130 {
2131 info = -2;
2132
2133 mattype.mark_as_unsymmetric ();
2134 typ = MatrixType::Full;
2135 }
2136 else
2137 {
2138 if (calc_cond)
2139 {
2140 Array<Complex> z (dim_vector (2 * nc, 1));
2141 Complex *pz = z.fortran_vec ();
2142 Array<double> rz (dim_vector (nc, 1));
2143 double *prz = rz.fortran_vec ();
2144
2145 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2146 nr, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2147 rcon, pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2148 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
2149
2150 if (info != 0)
2151 info = -2;
2152
2153 volatile double rcond_plus_one = rcon + 1.0;
2154
2155 if (rcond_plus_one == 1.0 || xisnan (rcon))
2156 {
2157 info = -2;
2158
2159 if (sing_handler)
2160 sing_handler (rcon);
2161 else
2162 (*current_liboctave_error_handler)
2163 ("matrix singular to machine precision, rcond = %g",
2164 rcon);
2165 }
2166 }
2167
2168 if (info == 0)
2169 {
2170 retval = b;
2171 Complex *result = retval.fortran_vec ();
2172
2173 octave_idx_type b_nc = b.cols ();
2174
2175 F77_XFCN (zpotrs, ZPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows
(), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2176 nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows
(), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2177 result, b.rows (), infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows
(), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2178 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows
(), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
2179 }
2180 else
2181 {
2182 mattype.mark_as_unsymmetric ();
2183 typ = MatrixType::Full;
2184 }
2185 }
2186 }
2187
2188 if (typ == MatrixType::Full)
2189 {
2190 info = 0;
2191
2192 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
2193 octave_idx_type *pipvt = ipvt.fortran_vec ();
2194
2195 ComplexMatrix atmp = *this;
2196 Complex *tmp_data = atmp.fortran_vec ();
2197
2198 Array<Complex> z (dim_vector (2 * nc, 1));
2199 Complex *pz = z.fortran_vec ();
2200 Array<double> rz (dim_vector (2 * nc, 1));
2201 double *prz = rz.fortran_vec ();
2202
2203 // Calculate the norm of the matrix, for later use.
2204 if (anorm < 0.)
2205 anorm = atmp.abs ().sum ().row (static_cast<octave_idx_type>(0))
2206 .max ();
2207
2208 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrf_ (nr, nr, tmp_data, nr, pipvt, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2209
2210 // Throw-away extra info LAPACK gives so as to not change output.
2211 rcon = 0.0;
2212 if (info != 0)
2213 {
2214 info = -2;
2215
2216 if (sing_handler)
2217 sing_handler (rcon);
2218 else
2219 (*current_liboctave_error_handler)
2220 ("matrix singular to machine precision");
2221
2222 mattype.mark_as_rectangular ();
2223 }
2224 else
2225 {
2226 if (calc_cond)
2227 {
2228 // Now calculate the condition number for
2229 // non-singular matrix.
2230 char job = '1';
2231 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2232 nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2233 rcon, pz, prz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2234 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, prz
, info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
2235
2236 if (info != 0)
2237 info = -2;
2238
2239 volatile double rcond_plus_one = rcon + 1.0;
2240
2241 if (rcond_plus_one == 1.0 || xisnan (rcon))
2242 {
2243 info = -2;
2244
2245 if (sing_handler)
2246 sing_handler (rcon);
2247 else
2248 (*current_liboctave_error_handler)
2249 ("matrix singular to machine precision, rcond = %g",
2250 rcon);
2251 }
2252 }
2253
2254 if (info == 0)
2255 {
2256 retval = b;
2257 Complex *result = retval.fortran_vec ();
2258
2259 octave_idx_type b_nc = b.cols ();
2260
2261 char job = 'N';
2262 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result,
b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2263 nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result,
b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2264 pipvt, result, b.rows (), infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result,
b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
2265 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result,
b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
2266 }
2267 else
2268 mattype.mark_as_rectangular ();
2269 }
2270 }
2271 }
2272
2273 return retval;
2274}
2275
2276ComplexMatrix
2277ComplexMatrix::solve (MatrixType &typ, const Matrix& b) const
2278{
2279 octave_idx_type info;
2280 double rcon;
2281 return solve (typ, b, info, rcon, 0);
2282}
2283
2284ComplexMatrix
2285ComplexMatrix::solve (MatrixType &typ, const Matrix& b,
2286 octave_idx_type& info) const
2287{
2288 double rcon;
2289 return solve (typ, b, info, rcon, 0);
2290}
2291
2292ComplexMatrix
2293ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
2294 double& rcon) const
2295{
2296 return solve (typ, b, info, rcon, 0);
2297}
2298
2299ComplexMatrix
2300ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
2301 double& rcon, solve_singularity_handler sing_handler,
2302 bool singular_fallback, blas_trans_type transt) const
2303{
2304 ComplexMatrix tmp (b);
2305 return solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt);
2306}
2307
2308ComplexMatrix
2309ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b) const
2310{
2311 octave_idx_type info;
2312 double rcon;
2313 return solve (typ, b, info, rcon, 0);
2314}
2315
2316ComplexMatrix
2317ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b,
2318 octave_idx_type& info) const
2319{
2320 double rcon;
2321 return solve (typ, b, info, rcon, 0);
2322}
2323
2324ComplexMatrix
2325ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b,
2326 octave_idx_type& info, double& rcon) const
2327{
2328 return solve (typ, b, info, rcon, 0);
2329}
2330
2331ComplexMatrix
2332ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
2333 octave_idx_type& info, double& rcon,
2334 solve_singularity_handler sing_handler,
2335 bool singular_fallback, blas_trans_type transt) const
2336{
2337 ComplexMatrix retval;
2338 int typ = mattype.type ();
2339
2340 if (typ == MatrixType::Unknown)
2341 typ = mattype.type (*this);
2342
2343 // Only calculate the condition number for LU/Cholesky
2344 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
2345 retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt);
2346 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
2347 retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt);
2348 else if (transt == blas_trans)
2349 return transpose ().solve (mattype, b, info, rcon, sing_handler,
2350 singular_fallback);
2351 else if (transt == blas_conj_trans)
2352 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler,
2353 singular_fallback);
2354 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
2355 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
2356 else if (typ != MatrixType::Rectangular)
2357 {
2358 (*current_liboctave_error_handler) ("unknown matrix type");
2359 return ComplexMatrix ();
2360 }
2361
2362 // Rectangular or one of the above solvers flags a singular matrix
2363 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
2364 {
2365 octave_idx_type rank;
2366 retval = lssolve (b, info, rank, rcon);
2367 }
2368
2369 return retval;
2370}
2371
2372ComplexColumnVector
2373ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b) const
2374{
2375 octave_idx_type info;
2376 double rcon;
2377 return solve (typ, ComplexColumnVector (b), info, rcon, 0);
2378}
2379
2380ComplexColumnVector
2381ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
2382 octave_idx_type& info) const
2383{
2384 double rcon;
2385 return solve (typ, ComplexColumnVector (b), info, rcon, 0);
2386}
2387
2388ComplexColumnVector
2389ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
2390 octave_idx_type& info, double& rcon) const
2391{
2392 return solve (typ, ComplexColumnVector (b), info, rcon, 0);
2393}
2394
2395ComplexColumnVector
2396ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
2397 octave_idx_type& info, double& rcon,
2398 solve_singularity_handler sing_handler,
2399 blas_trans_type transt) const
2400{
2401 return solve (typ, ComplexColumnVector (b), info, rcon, sing_handler, transt);
2402}
2403
2404ComplexColumnVector
2405ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b) const
2406{
2407 octave_idx_type info;
2408 double rcon;
2409 return solve (typ, b, info, rcon, 0);
2410}
2411
2412ComplexColumnVector
2413ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
2414 octave_idx_type& info) const
2415{
2416 double rcon;
2417 return solve (typ, b, info, rcon, 0);
2418}
2419
2420ComplexColumnVector
2421ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
2422 octave_idx_type& info, double& rcon) const
2423{
2424 return solve (typ, b, info, rcon, 0);
2425}
2426
2427ComplexColumnVector
2428ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
2429 octave_idx_type& info, double& rcon,
2430 solve_singularity_handler sing_handler,
2431 blas_trans_type transt) const
2432{
2433
2434 ComplexMatrix tmp (b);
2435 tmp = solve (typ, tmp, info, rcon, sing_handler, true, transt);
2436 return tmp.column (static_cast<octave_idx_type> (0));
2437}
2438
2439ComplexMatrix
2440ComplexMatrix::solve (const Matrix& b) const
2441{
2442 octave_idx_type info;
2443 double rcon;
2444 return solve (b, info, rcon, 0);
2445}
2446
2447ComplexMatrix
2448ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const
2449{
2450 double rcon;
2451 return solve (b, info, rcon, 0);
2452}
2453
2454ComplexMatrix
2455ComplexMatrix::solve (const Matrix& b, octave_idx_type& info,
2456 double& rcon) const
2457{
2458 return solve (b, info, rcon, 0);
2459}
2460
2461ComplexMatrix
2462ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon,
2463 solve_singularity_handler sing_handler,
2464 blas_trans_type transt) const
2465{
2466 ComplexMatrix tmp (b);
2467 return solve (tmp, info, rcon, sing_handler, transt);
2468}
2469
2470ComplexMatrix
2471ComplexMatrix::solve (const ComplexMatrix& b) const
2472{
2473 octave_idx_type info;
2474 double rcon;
2475 return solve (b, info, rcon, 0);
2476}
2477
2478ComplexMatrix
2479ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const
2480{
2481 double rcon;
2482 return solve (b, info, rcon, 0);
2483}
2484
2485ComplexMatrix
2486ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info,
2487 double& rcon) const
2488{
2489 return solve (b, info, rcon, 0);
2490}
2491
2492ComplexMatrix
2493ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info,
2494 double& rcon,
2495 solve_singularity_handler sing_handler,
2496 blas_trans_type transt) const
2497{
2498 MatrixType mattype (*this);
2499 return solve (mattype, b, info, rcon, sing_handler, true, transt);
2500}
2501
2502ComplexColumnVector
2503ComplexMatrix::solve (const ColumnVector& b) const
2504{
2505 octave_idx_type info;
2506 double rcon;
2507 return solve (ComplexColumnVector (b), info, rcon, 0);
2508}
2509
2510ComplexColumnVector
2511ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info) const
2512{
2513 double rcon;
2514 return solve (ComplexColumnVector (b), info, rcon, 0);
2515}
2516
2517ComplexColumnVector
2518ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
2519 double& rcon) const
2520{
2521 return solve (ComplexColumnVector (b), info, rcon, 0);
2522}
2523
2524ComplexColumnVector
2525ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
2526 double& rcon,
2527 solve_singularity_handler sing_handler,
2528 blas_trans_type transt) const
2529{
2530 return solve (ComplexColumnVector (b), info, rcon, sing_handler, transt);
2531}
2532
2533ComplexColumnVector
2534ComplexMatrix::solve (const ComplexColumnVector& b) const
2535{
2536 octave_idx_type info;
2537 double rcon;
2538 return solve (b, info, rcon, 0);
2539}
2540
2541ComplexColumnVector
2542ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const
2543{
2544 double rcon;
2545 return solve (b, info, rcon, 0);
2546}
2547
2548ComplexColumnVector
2549ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
2550 double& rcon) const
2551{
2552 return solve (b, info, rcon, 0);
2553}
2554
2555ComplexColumnVector
2556ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
2557 double& rcon,
2558 solve_singularity_handler sing_handler,
2559 blas_trans_type transt) const
2560{
2561 MatrixType mattype (*this);
2562 return solve (mattype, b, info, rcon, sing_handler, transt);
2563}
2564
2565ComplexMatrix
2566ComplexMatrix::lssolve (const Matrix& b) const
2567{
2568 octave_idx_type info;
2569 octave_idx_type rank;
2570 double rcon;
2571 return lssolve (ComplexMatrix (b), info, rank, rcon);
2572}
2573
2574ComplexMatrix
2575ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info) const
2576{
2577 octave_idx_type rank;
2578 double rcon;
2579 return lssolve (ComplexMatrix (b), info, rank, rcon);
2580}
2581
2582ComplexMatrix
2583ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info,
2584 octave_idx_type& rank) const
2585{
2586 double rcon;
2587 return lssolve (ComplexMatrix (b), info, rank, rcon);
2588}
2589
2590ComplexMatrix
2591ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info,
2592 octave_idx_type& rank, double& rcon) const
2593{
2594 return lssolve (ComplexMatrix (b), info, rank, rcon);
2595}
2596
2597ComplexMatrix
2598ComplexMatrix::lssolve (const ComplexMatrix& b) const
2599{
2600 octave_idx_type info;
2601 octave_idx_type rank;
2602 double rcon;
2603 return lssolve (b, info, rank, rcon);
2604}
2605
2606ComplexMatrix
2607ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const
2608{
2609 octave_idx_type rank;
2610 double rcon;
2611 return lssolve (b, info, rank, rcon);
2612}
2613
2614ComplexMatrix
2615ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
2616 octave_idx_type& rank) const
2617{
2618 double rcon;
2619 return lssolve (b, info, rank, rcon);
2620}
2621
2622ComplexMatrix
2623ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
2624 octave_idx_type& rank, double& rcon) const
2625{
2626 ComplexMatrix retval;
2627
2628 octave_idx_type nrhs = b.cols ();
2629
2630 octave_idx_type m = rows ();
2631 octave_idx_type n = cols ();
2632
2633 if (m != b.rows ())
2634 (*current_liboctave_error_handler)
2635 ("matrix dimension mismatch solution of linear equations");
2636 else if (m== 0 || n == 0 || b.cols () == 0)
2637 retval = ComplexMatrix (n, b.cols (), Complex (0.0, 0.0));
2638 else
2639 {
2640 volatile octave_idx_type minmn = (m < n ? m : n);
2641 octave_idx_type maxmn = m > n ? m : n;
2642 rcon = -1.0;
2643
2644 if (m != n)
2645 {
2646 retval = ComplexMatrix (maxmn, nrhs);
2647
2648 for (octave_idx_type j = 0; j < nrhs; j++)
2649 for (octave_idx_type i = 0; i < m; i++)
2650 retval.elem (i, j) = b.elem (i, j);
2651 }
2652 else
2653 retval = b;
2654
2655 ComplexMatrix atmp = *this;
2656 Complex *tmp_data = atmp.fortran_vec ();
2657
2658 Complex *pretval = retval.fortran_vec ();
2659 Array<double> s (dim_vector (minmn, 1));
2660 double *ps = s.fortran_vec ();
2661
2662 // Ask ZGELSD what the dimension of WORK should be.
2663 octave_idx_type lwork = -1;
2664
2665 Array<Complex> work (dim_vector (1, 1));
2666
2667 octave_idx_type smlsiz;
2668 F77_FUNC (xilaenv, XILAENV)xilaenv_ (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6)"ZGELSD",
2669 F77_CONST_CHAR_ARG2 (" ", 1)" ",
2670 0, 0, 0, 0, smlsiz
2671 F77_CHAR_ARG_LEN (6), 6
2672 F77_CHAR_ARG_LEN (1), 1);
2673
2674 octave_idx_type mnthr;
2675 F77_FUNC (xilaenv, XILAENV)xilaenv_ (6, F77_CONST_CHAR_ARG2 ("ZGELSD", 6)"ZGELSD",
2676 F77_CONST_CHAR_ARG2 (" ", 1)" ",
2677 m, n, nrhs, -1, mnthr
2678 F77_CHAR_ARG_LEN (6), 6
2679 F77_CHAR_ARG_LEN (1), 1);
2680
2681 // We compute the size of rwork and iwork because ZGELSD in
2682 // older versions of LAPACK does not return them on a query
2683 // call.
2684 double dminmn = static_cast<double> (minmn);
2685 double dsmlsizp1 = static_cast<double> (smlsiz+1);
2686#if defined (HAVE_LOG21)
2687 double tmp = log2 (dminmn / dsmlsizp1);
2688#else
2689 double tmp = log (dminmn / dsmlsizp1) / log (2.0);
2690#endif
2691 octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1;
2692 if (nlvl < 0)
2693 nlvl = 0;
2694
2695 octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
2696 + 3*smlsiz*nrhs
2697 + std::max ((smlsiz+1)*(smlsiz+1),
2698 n*(1+nrhs) + 2*nrhs);
2699 if (lrwork < 1)
2700 lrwork = 1;
2701 Array<double> rwork (dim_vector (lrwork, 1));
2702 double *prwork = rwork.fortran_vec ();
2703
2704 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2705 if (liwork < 1)
2706 liwork = 1;
2707 Array<octave_idx_type> iwork (dim_vector (liwork, 1));
2708 octave_idx_type* piwork = iwork.fortran_vec ();
2709
2710 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2711 ps, rcon, rank, work.fortran_vec (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2712 lwork, prwork, piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2713
2714 // The workspace query is broken in at least LAPACK 3.0.0
2715 // through 3.1.1 when n >= mnthr. The obtuse formula below
2716 // should provide sufficient workspace for ZGELSD to operate
2717 // efficiently.
2718 if (n > m && n >= mnthr)
2719 {
2720 octave_idx_type addend = m;
2721
2722 if (2*m-4 > addend)
2723 addend = 2*m-4;
2724
2725 if (nrhs > addend)
2726 addend = nrhs;
2727
2728 if (n-3*m > addend)
2729 addend = n-3*m;
2730
2731 const octave_idx_type lworkaround = 4*m + m*m + addend;
2732
2733 if (std::real (work(0)) < lworkaround)
2734 work(0) = lworkaround;
2735 }
2736 else if (m >= n)
2737 {
2738 octave_idx_type lworkaround = 2*m + m*nrhs;
2739
2740 if (std::real (work(0)) < lworkaround)
2741 work(0) = lworkaround;
2742 }
2743
2744 lwork = static_cast<octave_idx_type> (std::real (work(0)));
2745 work.resize (dim_vector (lwork, 1));
2746
2747 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2748 maxmn, ps, rcon, rank,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2749 work.fortran_vec (), lwork,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2750 prwork, piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2751
2752 if (s.elem (0) == 0.0)
2753 rcon = 0.0;
2754 else
2755 rcon = s.elem (minmn - 1) / s.elem (0);
2756
2757 retval.resize (n, nrhs);
2758 }
2759
2760 return retval;
2761}
2762
2763ComplexColumnVector
2764ComplexMatrix::lssolve (const ColumnVector& b) const
2765{
2766 octave_idx_type info;
2767 octave_idx_type rank;
2768 double rcon;
2769 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2770}
2771
2772ComplexColumnVector
2773ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info) const
2774{
2775 octave_idx_type rank;
2776 double rcon;
2777 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2778}
2779
2780ComplexColumnVector
2781ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info,
2782 octave_idx_type& rank) const
2783{
2784 double rcon;
2785 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2786}
2787
2788ComplexColumnVector
2789ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info,
2790 octave_idx_type& rank, double& rcon) const
2791{
2792 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2793}
2794
2795ComplexColumnVector
2796ComplexMatrix::lssolve (const ComplexColumnVector& b) const
2797{
2798 octave_idx_type info;
2799 octave_idx_type rank;
2800 double rcon;
2801 return lssolve (b, info, rank, rcon);
2802}
2803
2804ComplexColumnVector
2805ComplexMatrix::lssolve (const ComplexColumnVector& b,
2806 octave_idx_type& info) const
2807{
2808 octave_idx_type rank;
2809 double rcon;
2810 return lssolve (b, info, rank, rcon);
2811}
2812
2813ComplexColumnVector
2814ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
2815 octave_idx_type& rank) const
2816{
2817 double rcon;
2818 return lssolve (b, info, rank, rcon);
2819
2820}
2821
2822ComplexColumnVector
2823ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
2824 octave_idx_type& rank, double& rcon) const
2825{
2826 ComplexColumnVector retval;
2827
2828 octave_idx_type nrhs = 1;
2829
2830 octave_idx_type m = rows ();
2831 octave_idx_type n = cols ();
2832
2833 if (m != b.length ())
2834 (*current_liboctave_error_handler)
2835 ("matrix dimension mismatch solution of linear equations");
2836 else if (m == 0 || n == 0 || b.cols () == 0)
2837 retval = ComplexColumnVector (n, Complex (0.0, 0.0));
2838 else
2839 {
2840 volatile octave_idx_type minmn = (m < n ? m : n);
2841 octave_idx_type maxmn = m > n ? m : n;
2842 rcon = -1.0;
2843
2844 if (m != n)
2845 {
2846 retval = ComplexColumnVector (maxmn);
2847
2848 for (octave_idx_type i = 0; i < m; i++)
2849 retval.elem (i) = b.elem (i);
2850 }
2851 else
2852 retval = b;
2853
2854 ComplexMatrix atmp = *this;
2855 Complex *tmp_data = atmp.fortran_vec ();
2856
2857 Complex *pretval = retval.fortran_vec ();
2858 Array<double> s (dim_vector (minmn, 1));
2859 double *ps = s.fortran_vec ();
2860
2861 // Ask ZGELSD what the dimension of WORK should be.
2862 octave_idx_type lwork = -1;
2863
2864 Array<Complex> work (dim_vector (1, 1));
2865
2866 octave_idx_type smlsiz;
2867 F77_FUNC (xilaenv, XILAENV)xilaenv_ (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6)"ZGELSD",
2868 F77_CONST_CHAR_ARG2 (" ", 1)" ",
2869 0, 0, 0, 0, smlsiz
2870 F77_CHAR_ARG_LEN (6), 6
2871 F77_CHAR_ARG_LEN (1), 1);
2872
2873 // We compute the size of rwork and iwork because ZGELSD in
2874 // older versions of LAPACK does not return them on a query
2875 // call.
2876 double dminmn = static_cast<double> (minmn);
2877 double dsmlsizp1 = static_cast<double> (smlsiz+1);
2878#if defined (HAVE_LOG21)
2879 double tmp = log2 (dminmn / dsmlsizp1);
2880#else
2881 double tmp = log (dminmn / dsmlsizp1) / log (2.0);
2882#endif
2883 octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1;
2884 if (nlvl < 0)
2885 nlvl = 0;
2886
2887 octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
2888 + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
2889 if (lrwork < 1)
2890 lrwork = 1;
2891 Array<double> rwork (dim_vector (lrwork, 1));
2892 double *prwork = rwork.fortran_vec ();
2893
2894 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2895 if (liwork < 1)
2896 liwork = 1;
2897 Array<octave_idx_type> iwork (dim_vector (liwork, 1));
2898 octave_idx_type* piwork = iwork.fortran_vec ();
2899
2900 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2901 ps, rcon, rank, work.fortran_vec (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2902 lwork, prwork, piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2903
2904 lwork = static_cast<octave_idx_type> (std::real (work(0)));
2905 work.resize (dim_vector (lwork, 1));
2906 rwork.resize (dim_vector (static_cast<octave_idx_type> (rwork(0)), 1));
2907 iwork.resize (dim_vector (iwork(0), 1));
2908
2909 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2910 maxmn, ps, rcon, rank,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2911 work.fortran_vec (), lwork,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
2912 prwork, piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon
, rank, work.fortran_vec (), lwork, prwork, piwork, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
2913
2914 if (rank < minmn)
2915 {
2916 if (s.elem (0) == 0.0)
2917 rcon = 0.0;
2918 else
2919 rcon = s.elem (minmn - 1) / s.elem (0);
2920
2921 retval.resize (n, nrhs);
2922 }
2923 }
2924
2925 return retval;
2926}
2927
2928// column vector by row vector -> matrix operations
2929
2930ComplexMatrix
2931operator * (const ColumnVector& v, const ComplexRowVector& a)
2932{
2933 ComplexColumnVector tmp (v);
2934 return tmp * a;
2935}
2936
2937ComplexMatrix
2938operator * (const ComplexColumnVector& a, const RowVector& b)
2939{
2940 ComplexRowVector tmp (b);
2941 return a * tmp;
2942}
2943
2944ComplexMatrix
2945operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
2946{
2947 ComplexMatrix retval;
2948
2949 octave_idx_type len = v.length ();
2950
2951 if (len != 0)
2952 {
2953 octave_idx_type a_len = a.length ();
2954
2955 retval = ComplexMatrix (len, a_len);
2956 Complex *c = retval.fortran_vec ();
2957
2958 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data
(), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--;
octave_restore_current_context (saved_context); } } while (0
)
2959 F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data
(), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--;
octave_restore_current_context (saved_context); } } while (0
)
2960 len, a_len, 1, 1.0, v.data (), len,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data
(), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--;
octave_restore_current_context (saved_context); } } while (0
)
2961 a.data (), 1, 0.0, c, lendo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data
(), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--;
octave_restore_current_context (saved_context); } } while (0
)
2962 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data
(), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--;
octave_restore_current_context (saved_context); } } while (0
)
2963 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data
(), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--;
octave_restore_current_context (saved_context); } } while (0
)
;
2964 }
2965
2966 return retval;
2967}
2968
2969// matrix by diagonal matrix -> matrix operations
2970
2971ComplexMatrix&
2972ComplexMatrix::operator += (const DiagMatrix& a)
2973{
2974 octave_idx_type nr = rows ();
2975 octave_idx_type nc = cols ();
2976
2977 octave_idx_type a_nr = rows ();
2978 octave_idx_type a_nc = cols ();
2979
2980 if (nr != a_nr || nc != a_nc)
2981 {
2982 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2983 return *this;
2984 }
2985
2986 for (octave_idx_type i = 0; i < a.length (); i++)
2987 elem (i, i) += a.elem (i, i);
2988
2989 return *this;
2990}
2991
2992ComplexMatrix&
2993ComplexMatrix::operator -= (const DiagMatrix& a)
2994{
2995 octave_idx_type nr = rows ();
2996 octave_idx_type nc = cols ();
2997
2998 octave_idx_type a_nr = rows ();
2999 octave_idx_type a_nc = cols ();
3000
3001 if (nr != a_nr || nc != a_nc)
3002 {
3003 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
3004 return *this;
3005 }
3006
3007 for (octave_idx_type i = 0; i < a.length (); i++)
3008 elem (i, i) -= a.elem (i, i);
3009
3010 return *this;
3011}
3012
3013ComplexMatrix&
3014ComplexMatrix::operator += (const ComplexDiagMatrix& a)
3015{
3016 octave_idx_type nr = rows ();
3017 octave_idx_type nc = cols ();
3018
3019 octave_idx_type a_nr = rows ();
3020 octave_idx_type a_nc = cols ();
3021
3022 if (nr != a_nr || nc != a_nc)
3023 {
3024 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
3025 return *this;
3026 }
3027
3028 for (octave_idx_type i = 0; i < a.length (); i++)
3029 elem (i, i) += a.elem (i, i);
3030
3031 return *this;
3032}
3033
3034ComplexMatrix&
3035ComplexMatrix::operator -= (const ComplexDiagMatrix& a)
3036{
3037 octave_idx_type nr = rows ();
3038 octave_idx_type nc = cols ();
3039
3040 octave_idx_type a_nr = rows ();
3041 octave_idx_type a_nc = cols ();
3042
3043 if (nr != a_nr || nc != a_nc)
3044 {
3045 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
3046 return *this;
3047 }
3048
3049 for (octave_idx_type i = 0; i < a.length (); i++)
3050 elem (i, i) -= a.elem (i, i);
3051
3052 return *this;
3053}
3054
3055// matrix by matrix -> matrix operations
3056
3057ComplexMatrix&
3058ComplexMatrix::operator += (const Matrix& a)
3059{
3060 octave_idx_type nr = rows ();
3061 octave_idx_type nc = cols ();
3062
3063 octave_idx_type a_nr = a.rows ();
3064 octave_idx_type a_nc = a.cols ();
3065
3066 if (nr != a_nr || nc != a_nc)
3067 {
3068 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
3069 return *this;
3070 }
3071
3072 if (nr == 0 || nc == 0)
3073 return *this;
3074
3075 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
3076
3077 mx_inline_add2 (length (), d, a.data ());
3078 return *this;
3079}
3080
3081ComplexMatrix&
3082ComplexMatrix::operator -= (const Matrix& a)
3083{
3084 octave_idx_type nr = rows ();
3085 octave_idx_type nc = cols ();
3086
3087 octave_idx_type a_nr = a.rows ();
3088 octave_idx_type a_nc = a.cols ();
3089
3090 if (nr != a_nr || nc != a_nc)
3091 {
3092 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
3093 return *this;
3094 }
3095
3096 if (nr == 0 || nc == 0)
3097 return *this;
3098
3099 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
3100
3101 mx_inline_sub2 (length (), d, a.data ());
3102 return *this;
3103}
3104
3105// unary operations
3106
3107boolMatrix
3108ComplexMatrix::operator ! (void) const
3109{
3110 if (any_element_is_nan ())
3111 gripe_nan_to_logical_conversion ();
3112
3113 return do_mx_unary_op<bool, Complex> (*this, mx_inline_not);
3114}
3115
3116// other operations
3117
3118bool
3119ComplexMatrix::any_element_is_nan (void) const
3120{
3121 return do_mx_check<Complex> (*this, mx_inline_any_nan);
3122}
3123
3124bool
3125ComplexMatrix::any_element_is_inf_or_nan (void) const
3126{
3127 return ! do_mx_check<Complex> (*this, mx_inline_all_finite);
3128}
3129
3130// Return true if no elements have imaginary components.
3131
3132bool
3133ComplexMatrix::all_elements_are_real (void) const
3134{
3135 return do_mx_check<Complex> (*this, mx_inline_all_real);
3136}
3137
3138// Return nonzero if any element of CM has a non-integer real or
3139// imaginary part. Also extract the largest and smallest (real or
3140// imaginary) values and return them in MAX_VAL and MIN_VAL.
3141
3142bool
3143ComplexMatrix::all_integers (double& max_val, double& min_val) const
3144{
3145 octave_idx_type nr = rows ();
3146 octave_idx_type nc = cols ();
3147
3148 if (nr > 0 && nc > 0)
3149 {
3150 Complex val = elem (0, 0);
3151
3152 double r_val = std::real (val);
3153 double i_val = std::imag (val);
3154
3155 max_val = r_val;
3156 min_val = r_val;
3157
3158 if (i_val > max_val)
3159 max_val = i_val;
3160
3161 if (i_val < max_val)
3162 min_val = i_val;
3163 }
3164 else
3165 return false;
3166
3167 for (octave_idx_type j = 0; j < nc; j++)
3168 for (octave_idx_type i = 0; i < nr; i++)
3169 {
3170 Complex val = elem (i, j);
3171
3172 double r_val = std::real (val);
3173 double i_val = std::imag (val);
3174
3175 if (r_val > max_val)
3176 max_val = r_val;
3177
3178 if (i_val > max_val)
3179 max_val = i_val;
3180
3181 if (r_val < min_val)
3182 min_val = r_val;
3183
3184 if (i_val < min_val)
3185 min_val = i_val;
3186
3187 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
3188 return false;
3189 }
3190
3191 return true;
3192}
3193
3194bool
3195ComplexMatrix::too_large_for_float (void) const
3196{
3197 return test_any (xtoo_large_for_float);
3198}
3199
3200// FIXME: Do these really belong here? Maybe they should be in a base class?
3201
3202boolMatrix
3203ComplexMatrix::all (int dim) const
3204{
3205 return do_mx_red_op<bool, Complex> (*this, dim, mx_inline_all);
3206}
3207
3208boolMatrix
3209ComplexMatrix::any (int dim) const
3210{
3211 return do_mx_red_op<bool, Complex> (*this, dim, mx_inline_any);
3212}
3213
3214ComplexMatrix
3215ComplexMatrix::cumprod (int dim) const
3216{
3217 return do_mx_cum_op<Complex, Complex> (*this, dim, mx_inline_cumprod);
3218}
3219
3220ComplexMatrix
3221ComplexMatrix::cumsum (int dim) const
3222{
3223 return do_mx_cum_op<Complex, Complex> (*this, dim, mx_inline_cumsum);
3224}
3225
3226ComplexMatrix
3227ComplexMatrix::prod (int dim) const
3228{
3229 return do_mx_red_op<Complex, Complex> (*this, dim, mx_inline_prod);
3230}
3231
3232ComplexMatrix
3233ComplexMatrix::sum (int dim) const
3234{
3235 return do_mx_red_op<Complex, Complex> (*this, dim, mx_inline_sum);
3236}
3237
3238ComplexMatrix
3239ComplexMatrix::sumsq (int dim) const
3240{
3241 return do_mx_red_op<double, Complex> (*this, dim, mx_inline_sumsq);
3242}
3243
3244Matrix ComplexMatrix::abs (void) const
3245{
3246 return do_mx_unary_map<double, Complex, std::abs> (*this);
3247}
3248
3249ComplexMatrix
3250ComplexMatrix::diag (octave_idx_type k) const
3251{
3252 return MArray<Complex>::diag (k);
3253}
3254
3255ComplexDiagMatrix
3256ComplexMatrix::diag (octave_idx_type m, octave_idx_type n) const
3257{
3258 ComplexDiagMatrix retval;
3259
3260 octave_idx_type nr = rows ();
3261 octave_idx_type nc = cols ();
3262
3263 if (nr == 1 || nc == 1)
3264 retval = ComplexDiagMatrix (*this, m, n);
3265 else
3266 (*current_liboctave_error_handler)
3267 ("diag: expecting vector argument");
3268
3269 return retval;
3270}
3271
3272bool
3273ComplexMatrix::row_is_real_only (octave_idx_type i) const
3274{
3275 bool retval = true;
3276
3277 octave_idx_type nc = columns ();
3278
3279 for (octave_idx_type j = 0; j < nc; j++)
3280 {
3281 if (std::imag (elem (i, j)) != 0.0)
3282 {
3283 retval = false;
3284 break;
3285 }
3286 }
3287
3288 return retval;
3289}
3290
3291bool
3292ComplexMatrix::column_is_real_only (octave_idx_type j) const
3293{
3294 bool retval = true;
3295
3296 octave_idx_type nr = rows ();
3297
3298 for (octave_idx_type i = 0; i < nr; i++)
3299 {
3300 if (std::imag (elem (i, j)) != 0.0)
3301 {
3302 retval = false;
3303 break;
3304 }
3305 }
3306
3307 return retval;
3308}
3309
3310ComplexColumnVector
3311ComplexMatrix::row_min (void) const
3312{
3313 Array<octave_idx_type> dummy_idx;
3314 return row_min (dummy_idx);
3315}
3316
3317ComplexColumnVector
3318ComplexMatrix::row_min (Array<octave_idx_type>& idx_arg) const
3319{
3320 ComplexColumnVector result;
3321
3322 octave_idx_type nr = rows ();
3323 octave_idx_type nc = cols ();
3324
3325 if (nr > 0 && nc > 0)
3326 {
3327 result.resize (nr);
3328 idx_arg.resize (dim_vector (nr, 1));
3329
3330 for (octave_idx_type i = 0; i < nr; i++)
3331 {
3332 bool real_only = row_is_real_only (i);
3333
3334 octave_idx_type idx_j;
3335
3336 Complex tmp_min;
3337
3338 double abs_min = octave_NaN;
3339
3340 for (idx_j = 0; idx_j < nc; idx_j++)
3341 {
3342 tmp_min = elem (i, idx_j);
3343
3344 if (! xisnan (tmp_min))
3345 {
3346 abs_min = real_only ? std::real (tmp_min)
3347 : std::abs (tmp_min);
3348 break;
3349 }
3350 }
3351
3352 for (octave_idx_type j = idx_j+1; j < nc; j++)
3353 {
3354 Complex tmp = elem (i, j);
3355
3356 if (xisnan (tmp))
3357 continue;
3358
3359 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
3360
3361 if (abs_tmp < abs_min)
3362 {
3363 idx_j = j;
3364 tmp_min = tmp;
3365 abs_min = abs_tmp;
3366 }
3367 }
3368
3369 if (xisnan (tmp_min))
3370 {
3371 result.elem (i) = Complex_NaN_result;
3372 idx_arg.elem (i) = 0;
3373 }
3374 else
3375 {
3376 result.elem (i) = tmp_min;
3377 idx_arg.elem (i) = idx_j;
3378 }
3379 }
3380 }
3381
3382 return result;
3383}
3384
3385ComplexColumnVector
3386ComplexMatrix::row_max (void) const
3387{
3388 Array<octave_idx_type> dummy_idx;
3389 return row_max (dummy_idx);
3390}
3391
3392ComplexColumnVector
3393ComplexMatrix::row_max (Array<octave_idx_type>& idx_arg) const
3394{
3395 ComplexColumnVector result;
3396
3397 octave_idx_type nr = rows ();
3398 octave_idx_type nc = cols ();
3399
3400 if (nr > 0 && nc > 0)
3401 {
3402 result.resize (nr);
3403 idx_arg.resize (dim_vector (nr, 1));
3404
3405 for (octave_idx_type i = 0; i < nr; i++)
3406 {
3407 bool real_only = row_is_real_only (i);
3408
3409 octave_idx_type idx_j;
3410
3411 Complex tmp_max;
3412
3413 double abs_max = octave_NaN;
3414
3415 for (idx_j = 0; idx_j < nc; idx_j++)
3416 {
3417 tmp_max = elem (i, idx_j);
3418
3419 if (! xisnan (tmp_max))
3420 {
3421 abs_max = real_only ? std::real (tmp_max)
3422 : std::abs (tmp_max);
3423 break;
3424 }
3425 }
3426
3427 for (octave_idx_type j = idx_j+1; j < nc; j++)
3428 {
3429 Complex tmp = elem (i, j);
3430
3431 if (xisnan (tmp))
3432 continue;
3433
3434 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
3435
3436 if (abs_tmp > abs_max)
3437 {
3438 idx_j = j;
3439 tmp_max = tmp;
3440 abs_max = abs_tmp;
3441 }
3442 }
3443
3444 if (xisnan (tmp_max))
3445 {
3446 result.elem (i) = Complex_NaN_result;
3447 idx_arg.elem (i) = 0;
3448 }
3449 else
3450 {
3451 result.elem (i) = tmp_max;
3452 idx_arg.elem (i) = idx_j;
3453 }
3454 }
3455 }
3456
3457 return result;
3458}
3459
3460ComplexRowVector
3461ComplexMatrix::column_min (void) const
3462{
3463 Array<octave_idx_type> dummy_idx;
3464 return column_min (dummy_idx);
3465}
3466
3467ComplexRowVector
3468ComplexMatrix::column_min (Array<octave_idx_type>& idx_arg) const
3469{
3470 ComplexRowVector result;
3471
3472 octave_idx_type nr = rows ();
3473 octave_idx_type nc = cols ();
3474
3475 if (nr > 0 && nc > 0)
3476 {
3477 result.resize (nc);
3478 idx_arg.resize (dim_vector (1, nc));
3479
3480 for (octave_idx_type j = 0; j < nc; j++)
3481 {
3482 bool real_only = column_is_real_only (j);
3483
3484 octave_idx_type idx_i;
3485
3486 Complex tmp_min;
3487
3488 double abs_min = octave_NaN;
3489
3490 for (idx_i = 0; idx_i < nr; idx_i++)
3491 {
3492 tmp_min = elem (idx_i, j);
3493
3494 if (! xisnan (tmp_min))
3495 {
3496 abs_min = real_only ? std::real (tmp_min)
3497 : std::abs (tmp_min);
3498 break;
3499 }
3500 }
3501
3502 for (octave_idx_type i = idx_i+1; i < nr; i++)
3503 {
3504 Complex tmp = elem (i, j);
3505
3506 if (xisnan (tmp))
3507 continue;
3508
3509 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
3510
3511 if (abs_tmp < abs_min)
3512 {
3513 idx_i = i;
3514 tmp_min = tmp;
3515 abs_min = abs_tmp;
3516 }
3517 }
3518
3519 if (xisnan (tmp_min))
3520 {
3521 result.elem (j) = Complex_NaN_result;
3522 idx_arg.elem (j) = 0;
3523 }
3524 else
3525 {
3526 result.elem (j) = tmp_min;
3527 idx_arg.elem (j) = idx_i;
3528 }
3529 }
3530 }
3531
3532 return result;
3533}
3534
3535ComplexRowVector
3536ComplexMatrix::column_max (void) const
3537{
3538 Array<octave_idx_type> dummy_idx;
3539 return column_max (dummy_idx);
3540}
3541
3542ComplexRowVector
3543ComplexMatrix::column_max (Array<octave_idx_type>& idx_arg) const
3544{
3545 ComplexRowVector result;
3546
3547 octave_idx_type nr = rows ();
3548 octave_idx_type nc = cols ();
3549
3550 if (nr > 0 && nc > 0)
3551 {
3552 result.resize (nc);
3553 idx_arg.resize (dim_vector (1, nc));
3554
3555 for (octave_idx_type j = 0; j < nc; j++)
3556 {
3557 bool real_only = column_is_real_only (j);
3558
3559 octave_idx_type idx_i;
3560
3561 Complex tmp_max;
3562
3563 double abs_max = octave_NaN;
3564
3565 for (idx_i = 0; idx_i < nr; idx_i++)
3566 {
3567 tmp_max = elem (idx_i, j);
3568
3569 if (! xisnan (tmp_max))
3570 {
3571 abs_max = real_only ? std::real (tmp_max)
3572 : std::abs (tmp_max);
3573 break;
3574 }
3575 }
3576
3577 for (octave_idx_type i = idx_i+1; i < nr; i++)
3578 {
3579 Complex tmp = elem (i, j);
3580
3581 if (xisnan (tmp))
3582 continue;
3583
3584 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp);
3585
3586 if (abs_tmp > abs_max)
3587 {
3588 idx_i = i;
3589 tmp_max = tmp;
3590 abs_max = abs_tmp;
3591 }
3592 }
3593
3594 if (xisnan (tmp_max))
3595 {
3596 result.elem (j) = Complex_NaN_result;
3597 idx_arg.elem (j) = 0;
3598 }
3599 else
3600 {
3601 result.elem (j) = tmp_max;
3602 idx_arg.elem (j) = idx_i;
3603 }
3604 }
3605 }
3606
3607 return result;
3608}
3609
3610// i/o
3611
3612std::ostream&
3613operator << (std::ostream& os, const ComplexMatrix& a)
3614{
3615 for (octave_idx_type i = 0; i < a.rows (); i++)
3616 {
3617 for (octave_idx_type j = 0; j < a.cols (); j++)
3618 {
3619 os << " ";
3620 octave_write_complex (os, a.elem (i, j));
3621 }
3622 os << "\n";
3623 }
3624 return os;
3625}
3626
3627std::istream&
3628operator >> (std::istream& is, ComplexMatrix& a)
3629{
3630 octave_idx_type nr = a.rows ();
3631 octave_idx_type nc = a.cols ();
3632
3633 if (nr > 0 && nc > 0)
3634 {
3635 Complex tmp;
3636 for (octave_idx_type i = 0; i < nr; i++)
3637 for (octave_idx_type j = 0; j < nc; j++)
3638 {
3639 tmp = octave_read_value<Complex> (is);
3640 if (is)
3641 a.elem (i, j) = tmp;
3642 else
3643 goto done;
3644 }
3645 }
3646
3647done:
3648
3649 return is;
3650}
3651
3652ComplexMatrix
3653Givens (const Complex& x, const Complex& y)
3654{
3655 double cc;
3656 Complex cs, temp_r;
3657
3658 F77_FUNC (zlartg, ZLARTG)zlartg_ (x, y, cc, cs, temp_r);
3659
3660 ComplexMatrix g (2, 2);
3661
3662 g.elem (0, 0) = cc;
3663 g.elem (1, 1) = cc;
3664 g.elem (0, 1) = cs;
3665 g.elem (1, 0) = -conj (cs);
3666
3667 return g;
3668}
3669
3670ComplexMatrix
3671Sylvester (const ComplexMatrix& a, const ComplexMatrix& b,
3672 const ComplexMatrix& c)
3673{
3674 ComplexMatrix retval;
3675
3676 // FIXME: need to check that a, b, and c are all the same size.
3677
3678 // Compute Schur decompositions
3679
3680 ComplexSCHUR as (a, "U");
3681 ComplexSCHUR bs (b, "U");
3682
3683 // Transform c to new coordinates.
3684
3685 ComplexMatrix ua = as.unitary_matrix ();
3686 ComplexMatrix sch_a = as.schur_matrix ();
3687
3688 ComplexMatrix ub = bs.unitary_matrix ();
3689 ComplexMatrix sch_b = bs.schur_matrix ();
3690
3691 ComplexMatrix cx = ua.hermitian () * c * ub;
3692
3693 // Solve the sylvester equation, back-transform, and return the solution.
3694
3695 octave_idx_type a_nr = a.rows ();
3696 octave_idx_type b_nr = b.rows ();
3697
3698 double scale;
3699 octave_idx_type info;
3700
3701 Complex *pa = sch_a.fortran_vec ();
3702 Complex *pb = sch_b.fortran_vec ();
3703 Complex *px = cx.fortran_vec ();
3704
3705 F77_XFCN (ztrsyl, ZTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px,
a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3706 F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px,
a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3707 1, a_nr, b_nr, pa, a_nr, pb,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px,
a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3708 b_nr, px, a_nr, scale, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px,
a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3709 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px,
a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3710 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "ztrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; ztrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px,
a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
3711
3712 // FIXME: check info?
3713
3714 retval = -ua * cx * ub.hermitian ();
3715
3716 return retval;
3717}
3718
3719ComplexMatrix
3720operator * (const ComplexMatrix& m, const Matrix& a)
3721{
3722 if (m.columns () > std::min (m.rows (), a.columns ()) / 10)
3723 return ComplexMatrix (real (m) * a, imag (m) * a);
3724 else
3725 return m * ComplexMatrix (a);
3726}
3727
3728ComplexMatrix
3729operator * (const Matrix& m, const ComplexMatrix& a)
3730{
3731 if (a.rows () > std::min (m.rows (), a.columns ()) / 10)
3732 return ComplexMatrix (m * real (a), m * imag (a));
3733 else
3734 return ComplexMatrix (m) * a;
3735}
3736
3737/*
3738
3739## Simple Dot Product, Matrix-Vector, and Matrix-Matrix Unit tests
3740%!assert ([1+i 2+i 3+i] * [ 4+i ; 5+i ; 6+i], 29+21i, 1e-14)
3741%!assert ([1+i 2+i ; 3+i 4+i ] * [5+i ; 6+i], [15 + 14i ; 37 + 18i], 1e-14)
3742%!assert ([1+i 2+i ; 3+i 4+i ] * [5+i 6+i ; 7+i 8+i], [17 + 15i 20 + 17i; 41 + 19i 48 + 21i], 1e-14)
3743%!assert ([1 i]*[i 0]', -i);
3744
3745## Test some simple identities
3746%!shared M, cv, rv
3747%! M = randn (10,10) + i*rand (10,10);
3748%! cv = randn (10,1) + i*rand (10,1);
3749%! rv = randn (1,10) + i*rand (1,10);
3750%!assert ([M*cv,M*cv], M*[cv,cv], 1e-14)
3751%!assert ([M.'*cv,M.'*cv], M.'*[cv,cv], 1e-14)
3752%!assert ([M'*cv,M'*cv], M'*[cv,cv], 1e-14)
3753%!assert ([rv*M;rv*M], [rv;rv]*M, 1e-14)
3754%!assert ([rv*M.';rv*M.'], [rv;rv]*M.', 1e-14)
3755%!assert ([rv*M';rv*M'], [rv;rv]*M', 1e-14)
3756%!assert (2*rv*cv, [rv,rv]*[cv;cv], 1e-14)
3757
3758*/
3759
3760static inline char
3761get_blas_trans_arg (bool trans, bool conj)
3762{
3763 return trans ? (conj ? 'C' : 'T') : 'N';
3764}
3765
3766// the general GEMM operation
3767
3768ComplexMatrix
3769xgemm (const ComplexMatrix& a, const ComplexMatrix& b,
3770 blas_trans_type transa, blas_trans_type transb)
3771{
3772 ComplexMatrix retval;
3773
3774 bool tra = transa != blas_no_trans, trb = transb != blas_no_trans;
3775 bool cja = transa == blas_conj_trans, cjb = transb == blas_conj_trans;
3776
3777 octave_idx_type a_nr = tra ? a.cols () : a.rows ();
3778 octave_idx_type a_nc = tra ? a.rows () : a.cols ();
3779
3780 octave_idx_type b_nr = trb ? b.cols () : b.rows ();
3781 octave_idx_type b_nc = trb ? b.rows () : b.cols ();
3782
3783 if (a_nc != b_nr)
3784 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
3785 else
3786 {
3787 if (a_nr == 0 || a_nc == 0 || b_nc == 0)
3788 retval = ComplexMatrix (a_nr, b_nc, 0.0);
3789 else if (a.data () == b.data () && a_nr == b_nc && tra != trb)
3790 {
3791 octave_idx_type lda = a.rows ();
3792
3793 // FIXME: looking at the reference BLAS, it appears that it
3794 // should not be necessary to initialize the output matrix if
3795 // BETA is 0 in the call to ZHERK, but ATLAS appears to
3796 // use the result matrix before zeroing the elements.
3797
3798 retval = ComplexMatrix (a_nr, b_nc, 0.0);
3799 Complex *c = retval.fortran_vec ();
3800
3801 const char ctra = get_blas_trans_arg (tra, cja);
3802 if (cja || cjb)
3803 {
3804 F77_XFCN (zherk, ZHERK, (F77_CONST_CHAR_ARG2 ("U", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zherk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zherk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3805 F77_CONST_CHAR_ARG2 (&ctra, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zherk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zherk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3806 a_nr, a_nc, 1.0,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zherk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zherk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3807 a.data (), lda, 0.0, c, a_nrdo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zherk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zherk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3808 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zherk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zherk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3809 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zherk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zherk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
3810 for (octave_idx_type j = 0; j < a_nr; j++)
3811 for (octave_idx_type i = 0; i < j; i++)
3812 retval.xelem (j,i) = std::conj (retval.xelem (i,j));
3813 }
3814 else
3815 {
3816 F77_XFCN (zsyrk, ZSYRK, (F77_CONST_CHAR_ARG2 ("U", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3817 F77_CONST_CHAR_ARG2 (&ctra, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3818 a_nr, a_nc, 1.0,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3819 a.data (), lda, 0.0, c, a_nrdo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3820 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3821 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0
, c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
3822 for (octave_idx_type j = 0; j < a_nr; j++)
3823 for (octave_idx_type i = 0; i < j; i++)
3824 retval.xelem (j,i) = retval.xelem (i,j);
3825
3826 }
3827
3828 }
3829 else
3830 {
3831 octave_idx_type lda = a.rows (), tda = a.cols ();
3832 octave_idx_type ldb = b.rows (), tdb = b.cols ();
3833
3834 retval = ComplexMatrix (a_nr, b_nc, 0.0);
3835 Complex *c = retval.fortran_vec ();
3836
3837 if (b_nc == 1 && a_nr == 1)
3838 {
3839 if (cja == cjb)
3840 {
3841 F77_FUNC (xzdotu, XZDOTU)xzdotu_ (a_nc, a.data (), 1, b.data (), 1,
3842 *c);
3843 if (cja) *c = std::conj (*c);
3844 }
3845 else if (cja)
3846 F77_FUNC (xzdotc, XZDOTC)xzdotc_ (a_nc, a.data (), 1, b.data (), 1,
3847 *c);
3848 else
3849 F77_FUNC (xzdotc, XZDOTC)xzdotc_ (a_nc, b.data (), 1, a.data (), 1,
3850 *c);
3851 }
3852 else if (b_nc == 1 && ! cjb)
3853 {
3854 const char ctra = get_blas_trans_arg (tra, cja);
3855 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 (&ctra, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data (
), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3856 lda, tda, 1.0, a.data (), lda,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data (
), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3857 b.data (), 1, 0.0, c, 1do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data (
), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3858 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data (
), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
3859 }
3860 else if (a_nr == 1 && ! cja && ! cjb)
3861 {
3862 const char crevtrb = get_blas_trans_arg (! trb, cjb);
3863 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 (&crevtrb, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data
(), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3864 ldb, tdb, 1.0, b.data (), ldb,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data
(), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3865 a.data (), 1, 0.0, c, 1do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data
(), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
3866 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data
(), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
3867 }
3868 else
3869 {
3870 const char ctra = get_blas_trans_arg (tra, cja);
3871 const char ctrb = get_blas_trans_arg (trb, cjb);
3872 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 (&ctra, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data
(), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
3873 F77_CONST_CHAR_ARG2 (&ctrb, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data
(), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
3874 a_nr, b_nc, a_nc, 1.0, a.data (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data
(), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
3875 lda, b.data (), ldb, 0.0, c, a_nrdo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data
(), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
3876 F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data
(), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
3877 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data
(), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
3878 }
3879 }
3880 }
3881
3882 return retval;
3883}
3884
3885ComplexMatrix
3886operator * (const ComplexMatrix& a, const ComplexMatrix& b)
3887{
3888 return xgemm (a, b);
3889}
3890
3891// FIXME: it would be nice to share code among the min/max functions below.
3892
3893#define EMPTY_RETURN_CHECK(T)if (nr == 0 || nc == 0) return T (nr, nc); \
3894 if (nr == 0 || nc == 0) \
3895 return T (nr, nc);
3896
3897ComplexMatrix
3898min (const Complex& c, const ComplexMatrix& m)
3899{
3900 octave_idx_type nr = m.rows ();
3901 octave_idx_type nc = m.columns ();
3902
3903 EMPTY_RETURN_CHECK (ComplexMatrix)if (nr == 0 || nc == 0) return ComplexMatrix (nr, nc);;
3904
3905 ComplexMatrix result (nr, nc);
3906
3907 for (octave_idx_type j = 0; j < nc; j++)
3908 for (octave_idx_type i = 0; i < nr; i++)
3909 {
3910 octave_quit ();
3911 result (i, j) = xmin (c, m (i, j));
3912 }
3913
3914 return result;
3915}
3916
3917ComplexMatrix
3918min (const ComplexMatrix& m, const Complex& c)
3919{
3920 octave_idx_type nr = m.rows ();
3921 octave_idx_type nc = m.columns ();
3922
3923 EMPTY_RETURN_CHECK (ComplexMatrix)if (nr == 0 || nc == 0) return ComplexMatrix (nr, nc);;
3924
3925 ComplexMatrix result (nr, nc);
3926
3927 for (octave_idx_type j = 0; j < nc; j++)
3928 for (octave_idx_type i = 0; i < nr; i++)
3929 {
3930 octave_quit ();
3931 result (i, j) = xmin (m (i, j), c);
3932 }
3933
3934 return result;
3935}
3936
3937ComplexMatrix
3938min (const ComplexMatrix& a, const ComplexMatrix& b)
3939{
3940 octave_idx_type nr = a.rows ();
3941 octave_idx_type nc = a.columns ();
3942
3943 if (nr != b.rows () || nc != b.columns ())
3944 {
3945 (*current_liboctave_error_handler)
3946 ("two-arg min expecting args of same size");
3947 return ComplexMatrix ();
3948 }
3949
3950 EMPTY_RETURN_CHECK (ComplexMatrix)if (nr == 0 || nc == 0) return ComplexMatrix (nr, nc);;
3951
3952 ComplexMatrix result (nr, nc);
3953
3954 for (octave_idx_type j = 0; j < nc; j++)
3955 {
3956 int columns_are_real_only = 1;
3957 for (octave_idx_type i = 0; i < nr; i++)
3958 {
3959 octave_quit ();
3960 if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0)
3961 {
3962 columns_are_real_only = 0;
3963 break;
3964 }
3965 }
3966
3967 if (columns_are_real_only)
3968 {
3969 for (octave_idx_type i = 0; i < nr; i++)
3970 result (i, j) = xmin (std::real (a (i, j)), std::real (b (i, j)));
3971 }
3972 else
3973 {
3974 for (octave_idx_type i = 0; i < nr; i++)
3975 {
3976 octave_quit ();
3977 result (i, j) = xmin (a (i, j), b (i, j));
3978 }
3979 }
3980 }
3981
3982 return result;
3983}
3984
3985ComplexMatrix
3986max (const Complex& c, const ComplexMatrix& m)
3987{
3988 octave_idx_type nr = m.rows ();
3989 octave_idx_type nc = m.columns ();
3990
3991 EMPTY_RETURN_CHECK (ComplexMatrix)if (nr == 0 || nc == 0) return ComplexMatrix (nr, nc);;
3992
3993 ComplexMatrix result (nr, nc);
3994
3995 for (octave_idx_type j = 0; j < nc; j++)
3996 for (octave_idx_type i = 0; i < nr; i++)
3997 {
3998 octave_quit ();
3999 result (i, j) = xmax (c, m (i, j));
4000 }
4001
4002 return result;
4003}
4004
4005ComplexMatrix
4006max (const ComplexMatrix& m, const Complex& c)
4007{
4008 octave_idx_type nr = m.rows ();
4009 octave_idx_type nc = m.columns ();
4010
4011 EMPTY_RETURN_CHECK (ComplexMatrix)if (nr == 0 || nc == 0) return ComplexMatrix (nr, nc);;
4012
4013 ComplexMatrix result (nr, nc);
4014
4015 for (octave_idx_type j = 0; j < nc; j++)
4016 for (octave_idx_type i = 0; i < nr; i++)
4017 {
4018 octave_quit ();
4019 result (i, j) = xmax (m (i, j), c);
4020 }
4021
4022 return result;
4023}
4024
4025ComplexMatrix
4026max (const ComplexMatrix& a, const ComplexMatrix& b)
4027{
4028 octave_idx_type nr = a.rows ();
4029 octave_idx_type nc = a.columns ();
4030
4031 if (nr != b.rows () || nc != b.columns ())
4032 {
4033 (*current_liboctave_error_handler)
4034 ("two-arg max expecting args of same size");
4035 return ComplexMatrix ();
4036 }
4037
4038 EMPTY_RETURN_CHECK (ComplexMatrix)if (nr == 0 || nc == 0) return ComplexMatrix (nr, nc);;
4039
4040 ComplexMatrix result (nr, nc);
4041
4042 for (octave_idx_type j = 0; j < nc; j++)
4043 {
4044 int columns_are_real_only = 1;
4045 for (octave_idx_type i = 0; i < nr; i++)
4046 {
4047 octave_quit ();
4048 if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0)
4049 {
4050 columns_are_real_only = 0;
4051 break;
4052 }
4053 }
4054
4055 if (columns_are_real_only)
4056 {
4057 for (octave_idx_type i = 0; i < nr; i++)
4058 {
4059 octave_quit ();
4060 result (i, j) = xmax (std::real (a (i, j)), std::real (b (i, j)));
4061 }
4062 }
4063 else
4064 {
4065 for (octave_idx_type i = 0; i < nr; i++)
4066 {
4067 octave_quit ();
4068 result (i, j) = xmax (a (i, j), b (i, j));
4069 }
4070 }
4071 }
4072
4073 return result;
4074}
4075
4076ComplexMatrix linspace (const ComplexColumnVector& x1,
4077 const ComplexColumnVector& x2,
4078 octave_idx_type n)
4079
4080{
4081 if (n < 1) n = 1;
4082
4083 octave_idx_type m = x1.length ();
4084
4085 if (x2.length () != m)
4086 (*current_liboctave_error_handler)
4087 ("linspace: vectors must be of equal length");
4088
4089 NoAlias<ComplexMatrix> retval;
4090
4091 retval.clear (m, n);
4092 for (octave_idx_type i = 0; i < m; i++)
4093 retval(i, 0) = x1(i);
4094
4095 // The last column is not needed while using delta.
4096 Complex *delta = &retval(0, n-1);
4097 for (octave_idx_type i = 0; i < m; i++)
4098 delta[i] = (x2(i) - x1(i)) / (n - 1.0);
4099
4100 for (octave_idx_type j = 1; j < n-1; j++)
4101 for (octave_idx_type i = 0; i < m; i++)
4102 retval(i, j) = x1(i) + static_cast<double> (j)*delta[i];
4103
4104 for (octave_idx_type i = 0; i < m; i++)
4105 retval(i, n-1) = x2(i);
4106
4107 return retval;
4108}
4109
4110MS_CMP_OPS (ComplexMatrix, Complex)boolMatrix mx_el_lt (const ComplexMatrix& m, const Complex
& s) { return do_ms_binary_op<bool, ComplexMatrix::element_type
, Complex> (m, s, mx_inline_lt); } boolMatrix mx_el_le (const
ComplexMatrix& m, const Complex& s) { return do_ms_binary_op
<bool, ComplexMatrix::element_type, Complex> (m, s, mx_inline_le
); } boolMatrix mx_el_ge (const ComplexMatrix& m, const Complex
& s) { return do_ms_binary_op<bool, ComplexMatrix::element_type
, Complex> (m, s, mx_inline_ge); } boolMatrix mx_el_gt (const
ComplexMatrix& m, const Complex& s) { return do_ms_binary_op
<bool, ComplexMatrix::element_type, Complex> (m, s, mx_inline_gt
); } boolMatrix mx_el_eq (const ComplexMatrix& m, const Complex
& s) { return do_ms_binary_op<bool, ComplexMatrix::element_type
, Complex> (m, s, mx_inline_eq); } boolMatrix mx_el_ne (const
ComplexMatrix& m, const Complex& s) { return do_ms_binary_op
<bool, ComplexMatrix::element_type, Complex> (m, s, mx_inline_ne
); }
4111MS_BOOL_OPS (ComplexMatrix, Complex)boolMatrix mx_el_and (const ComplexMatrix& m, const Complex
& s) { if (do_mx_check (m, mx_inline_any_nan<ComplexMatrix
::element_type>)) gripe_nan_to_logical_conversion (); if (
xisnan (s)) gripe_nan_to_logical_conversion (); return do_ms_binary_op
<bool, ComplexMatrix::element_type, Complex> (m, s, mx_inline_and
); } boolMatrix mx_el_or (const ComplexMatrix& m, const Complex
& s) { if (do_mx_check (m, mx_inline_any_nan<ComplexMatrix
::element_type>)) gripe_nan_to_logical_conversion (); if (
xisnan (s)) gripe_nan_to_logical_conversion (); return do_ms_binary_op
<bool, ComplexMatrix::element_type, Complex> (m, s, mx_inline_or
); }
4112
4113SM_CMP_OPS (Complex, ComplexMatrix)boolMatrix mx_el_lt (const Complex& s, const ComplexMatrix
& m) { return do_sm_binary_op<bool, Complex, ComplexMatrix
::element_type> (s, m, mx_inline_lt); } boolMatrix mx_el_le
(const Complex& s, const ComplexMatrix& m) { return do_sm_binary_op
<bool, Complex, ComplexMatrix::element_type> (s, m, mx_inline_le
); } boolMatrix mx_el_ge (const Complex& s, const ComplexMatrix
& m) { return do_sm_binary_op<bool, Complex, ComplexMatrix
::element_type> (s, m, mx_inline_ge); } boolMatrix mx_el_gt
(const Complex& s, const ComplexMatrix& m) { return do_sm_binary_op
<bool, Complex, ComplexMatrix::element_type> (s, m, mx_inline_gt
); } boolMatrix mx_el_eq (const Complex& s, const ComplexMatrix
& m) { return do_sm_binary_op<bool, Complex, ComplexMatrix
::element_type> (s, m, mx_inline_eq); } boolMatrix mx_el_ne
(const Complex& s, const ComplexMatrix& m) { return do_sm_binary_op
<bool, Complex, ComplexMatrix::element_type> (s, m, mx_inline_ne
); }
4114SM_BOOL_OPS (Complex, ComplexMatrix)boolMatrix mx_el_and (const Complex& s, const ComplexMatrix
& m) { if (xisnan (s)) gripe_nan_to_logical_conversion ()
; if (do_mx_check (m, mx_inline_any_nan<ComplexMatrix::element_type
>)) gripe_nan_to_logical_conversion (); return do_sm_binary_op
<bool, Complex, ComplexMatrix::element_type> (s, m, mx_inline_and
); } boolMatrix mx_el_or (const Complex& s, const ComplexMatrix
& m) { if (xisnan (s)) gripe_nan_to_logical_conversion ()
; if (do_mx_check (m, mx_inline_any_nan<ComplexMatrix::element_type
>)) gripe_nan_to_logical_conversion (); return do_sm_binary_op
<bool, Complex, ComplexMatrix::element_type> (s, m, mx_inline_or
); }
4115
4116MM_CMP_OPS (ComplexMatrix, ComplexMatrix)boolMatrix mx_el_lt (const ComplexMatrix& m1, const ComplexMatrix
& m2) { return do_mm_binary_op<bool, ComplexMatrix::element_type
, ComplexMatrix::element_type> (m1, m2, mx_inline_lt, mx_inline_lt
, mx_inline_lt, "mx_el_lt"); } boolMatrix mx_el_le (const ComplexMatrix
& m1, const ComplexMatrix& m2) { return do_mm_binary_op
<bool, ComplexMatrix::element_type, ComplexMatrix::element_type
> (m1, m2, mx_inline_le, mx_inline_le, mx_inline_le, "mx_el_le"
); } boolMatrix mx_el_ge (const ComplexMatrix& m1, const ComplexMatrix
& m2) { return do_mm_binary_op<bool, ComplexMatrix::element_type
, ComplexMatrix::element_type> (m1, m2, mx_inline_ge, mx_inline_ge
, mx_inline_ge, "mx_el_ge"); } boolMatrix mx_el_gt (const ComplexMatrix
& m1, const ComplexMatrix& m2) { return do_mm_binary_op
<bool, ComplexMatrix::element_type, ComplexMatrix::element_type
> (m1, m2, mx_inline_gt, mx_inline_gt, mx_inline_gt, "mx_el_gt"
); } boolMatrix mx_el_eq (const ComplexMatrix& m1, const ComplexMatrix
& m2) { return do_mm_binary_op<bool, ComplexMatrix::element_type
, ComplexMatrix::element_type> (m1, m2, mx_inline_eq, mx_inline_eq
, mx_inline_eq, "mx_el_eq"); } boolMatrix mx_el_ne (const ComplexMatrix
& m1, const ComplexMatrix& m2) { return do_mm_binary_op
<bool, ComplexMatrix::element_type, ComplexMatrix::element_type
> (m1, m2, mx_inline_ne, mx_inline_ne, mx_inline_ne, "mx_el_ne"
); }
4117MM_BOOL_OPS (ComplexMatrix, ComplexMatrix)boolMatrix mx_el_and (const ComplexMatrix& m1, const ComplexMatrix
& m2) { if (do_mx_check (m1, mx_inline_any_nan<ComplexMatrix
::element_type>)) gripe_nan_to_logical_conversion (); if (
do_mx_check (m2, mx_inline_any_nan<ComplexMatrix::element_type
>)) gripe_nan_to_logical_conversion (); return do_mm_binary_op
<bool, ComplexMatrix::element_type, ComplexMatrix::element_type
> (m1, m2, mx_inline_and, mx_inline_and, mx_inline_and, "mx_el_and"
); } boolMatrix mx_el_or (const ComplexMatrix& m1, const ComplexMatrix
& m2) { if (do_mx_check (m1, mx_inline_any_nan<ComplexMatrix
::element_type>)) gripe_nan_to_logical_conversion (); if (
do_mx_check (m2, mx_inline_any_nan<ComplexMatrix::element_type
>)) gripe_nan_to_logical_conversion (); return do_mm_binary_op
<bool, ComplexMatrix::element_type, ComplexMatrix::element_type
> (m1, m2, mx_inline_or, mx_inline_or, mx_inline_or, "mx_el_or"
); }