Bug Summary

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