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