|
1
|
1 /* |
|
|
2 |
|
2847
|
3 Copyright (C) 1996, 1997 John W. Eaton |
|
1
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
|
10 later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
19 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
20 02110-1301, USA. |
|
1
|
21 |
|
|
22 */ |
|
|
23 |
|
240
|
24 #ifdef HAVE_CONFIG_H |
|
1192
|
25 #include <config.h> |
|
1
|
26 #endif |
|
|
27 |
|
1343
|
28 #include <cassert> |
|
|
29 |
|
4669
|
30 #include "Array-util.h" |
|
1352
|
31 #include "CMatrix.h" |
|
453
|
32 #include "dMatrix.h" |
|
4543
|
33 #include "CNDArray.h" |
|
|
34 #include "dNDArray.h" |
|
1651
|
35 #include "oct-cmplx.h" |
|
4153
|
36 #include "quit.h" |
|
1352
|
37 |
|
|
38 #include "error.h" |
|
|
39 #include "xdiv.h" |
|
1
|
40 |
|
3480
|
41 static inline bool |
|
5275
|
42 result_ok (octave_idx_type info) |
|
1
|
43 { |
|
|
44 assert (info != -1); |
|
|
45 |
|
3480
|
46 return (info != -2); |
|
|
47 } |
|
932
|
48 |
|
3480
|
49 static void |
|
|
50 solve_singularity_warning (double rcond) |
|
|
51 { |
|
|
52 warning ("matrix singular to machine precision, rcond = %g", rcond); |
|
3719
|
53 warning ("attempting to find minimum norm solution"); |
|
1
|
54 } |
|
|
55 |
|
2364
|
56 template <class T1, class T2> |
|
|
57 bool |
|
3195
|
58 mx_leftdiv_conform (const T1& a, const T2& b) |
|
1
|
59 { |
|
5275
|
60 octave_idx_type a_nr = a.rows (); |
|
|
61 octave_idx_type b_nr = b.rows (); |
|
2364
|
62 |
|
1
|
63 if (a_nr != b_nr) |
|
|
64 { |
|
5275
|
65 octave_idx_type a_nc = a.cols (); |
|
|
66 octave_idx_type b_nc = b.cols (); |
|
2364
|
67 |
|
|
68 gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc); |
|
|
69 return false; |
|
1
|
70 } |
|
|
71 |
|
2364
|
72 return true; |
|
1
|
73 } |
|
|
74 |
|
3195
|
75 #define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \ |
|
|
76 template bool mx_leftdiv_conform (const T1&, const T2&) |
|
|
77 |
|
|
78 INSTANTIATE_MX_LEFTDIV_CONFORM (Matrix, Matrix); |
|
|
79 INSTANTIATE_MX_LEFTDIV_CONFORM (Matrix, ComplexMatrix); |
|
|
80 INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexMatrix, Matrix); |
|
|
81 INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexMatrix, ComplexMatrix); |
|
2364
|
82 |
|
|
83 template <class T1, class T2> |
|
|
84 bool |
|
3195
|
85 mx_div_conform (const T1& a, const T2& b) |
|
1
|
86 { |
|
5275
|
87 octave_idx_type a_nc = a.cols (); |
|
|
88 octave_idx_type b_nc = b.cols (); |
|
2364
|
89 |
|
1
|
90 if (a_nc != b_nc) |
|
|
91 { |
|
5275
|
92 octave_idx_type a_nr = a.rows (); |
|
|
93 octave_idx_type b_nr = b.rows (); |
|
2364
|
94 |
|
|
95 gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc); |
|
|
96 return false; |
|
1
|
97 } |
|
|
98 |
|
2364
|
99 return true; |
|
1
|
100 } |
|
|
101 |
|
3195
|
102 #define INSTANTIATE_MX_DIV_CONFORM(T1, T2) \ |
|
|
103 template bool mx_div_conform (const T1&, const T2&) |
|
|
104 |
|
|
105 INSTANTIATE_MX_DIV_CONFORM (Matrix, Matrix); |
|
|
106 INSTANTIATE_MX_DIV_CONFORM (Matrix, ComplexMatrix); |
|
|
107 INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, Matrix); |
|
|
108 INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, ComplexMatrix); |
|
2364
|
109 |
|
767
|
110 // Right division functions. |
|
|
111 // |
|
|
112 // op2 / op1: m cm |
|
|
113 // +-- +---+----+ |
|
|
114 // matrix | 1 | 3 | |
|
|
115 // +---+----+ |
|
|
116 // complex_matrix | 2 | 4 | |
|
|
117 // +---+----+ |
|
1
|
118 |
|
767
|
119 // -*- 1 -*- |
|
1800
|
120 Matrix |
|
164
|
121 xdiv (const Matrix& a, const Matrix& b) |
|
1
|
122 { |
|
2364
|
123 if (! mx_div_conform (a, b)) |
|
1800
|
124 return Matrix (); |
|
1
|
125 |
|
|
126 Matrix atmp = a.transpose (); |
|
|
127 Matrix btmp = b.transpose (); |
|
|
128 |
|
5275
|
129 octave_idx_type info; |
|
1
|
130 if (btmp.rows () == btmp.columns ()) |
|
|
131 { |
|
|
132 double rcond = 0.0; |
|
3480
|
133 |
|
|
134 Matrix result |
|
|
135 = btmp.solve (atmp, info, rcond, solve_singularity_warning); |
|
|
136 |
|
|
137 if (result_ok (info)) |
|
1800
|
138 return Matrix (result.transpose ()); |
|
1
|
139 } |
|
|
140 |
|
5275
|
141 octave_idx_type rank; |
|
1
|
142 Matrix result = btmp.lssolve (atmp, info, rank); |
|
|
143 |
|
1800
|
144 return result.transpose (); |
|
1
|
145 } |
|
|
146 |
|
767
|
147 // -*- 2 -*- |
|
1800
|
148 ComplexMatrix |
|
164
|
149 xdiv (const Matrix& a, const ComplexMatrix& b) |
|
1
|
150 { |
|
2364
|
151 if (! mx_div_conform (a, b)) |
|
1800
|
152 return ComplexMatrix (); |
|
1
|
153 |
|
|
154 Matrix atmp = a.transpose (); |
|
|
155 ComplexMatrix btmp = b.hermitian (); |
|
|
156 |
|
5275
|
157 octave_idx_type info; |
|
1
|
158 if (btmp.rows () == btmp.columns ()) |
|
|
159 { |
|
|
160 double rcond = 0.0; |
|
3480
|
161 |
|
|
162 ComplexMatrix result |
|
|
163 = btmp.solve (atmp, info, rcond, solve_singularity_warning); |
|
|
164 |
|
|
165 if (result_ok (info)) |
|
1800
|
166 return result.hermitian (); |
|
1
|
167 } |
|
|
168 |
|
5275
|
169 octave_idx_type rank; |
|
1
|
170 ComplexMatrix result = btmp.lssolve (atmp, info, rank); |
|
|
171 |
|
1800
|
172 return result.hermitian (); |
|
1
|
173 } |
|
|
174 |
|
767
|
175 // -*- 3 -*- |
|
1800
|
176 ComplexMatrix |
|
164
|
177 xdiv (const ComplexMatrix& a, const Matrix& b) |
|
1
|
178 { |
|
2364
|
179 if (! mx_div_conform (a, b)) |
|
1800
|
180 return ComplexMatrix (); |
|
1
|
181 |
|
|
182 ComplexMatrix atmp = a.hermitian (); |
|
|
183 Matrix btmp = b.transpose (); |
|
|
184 |
|
5275
|
185 octave_idx_type info; |
|
1
|
186 if (btmp.rows () == btmp.columns ()) |
|
|
187 { |
|
|
188 double rcond = 0.0; |
|
3480
|
189 |
|
|
190 ComplexMatrix result |
|
|
191 = btmp.solve (atmp, info, rcond, solve_singularity_warning); |
|
|
192 |
|
|
193 if (result_ok (info)) |
|
1800
|
194 return result.hermitian (); |
|
1
|
195 } |
|
|
196 |
|
5275
|
197 octave_idx_type rank; |
|
1
|
198 ComplexMatrix result = btmp.lssolve (atmp, info, rank); |
|
|
199 |
|
1800
|
200 return result.hermitian (); |
|
1
|
201 } |
|
|
202 |
|
767
|
203 // -*- 4 -*- |
|
1800
|
204 ComplexMatrix |
|
164
|
205 xdiv (const ComplexMatrix& a, const ComplexMatrix& b) |
|
1
|
206 { |
|
2364
|
207 if (! mx_div_conform (a, b)) |
|
1800
|
208 return ComplexMatrix (); |
|
1
|
209 |
|
|
210 ComplexMatrix atmp = a.hermitian (); |
|
|
211 ComplexMatrix btmp = b.hermitian (); |
|
|
212 |
|
5275
|
213 octave_idx_type info; |
|
1
|
214 if (btmp.rows () == btmp.columns ()) |
|
|
215 { |
|
|
216 double rcond = 0.0; |
|
3480
|
217 |
|
|
218 ComplexMatrix result |
|
|
219 = btmp.solve (atmp, info, rcond, solve_singularity_warning); |
|
|
220 |
|
|
221 if (result_ok (info)) |
|
1800
|
222 return result.hermitian (); |
|
1
|
223 } |
|
|
224 |
|
5275
|
225 octave_idx_type rank; |
|
1
|
226 ComplexMatrix result = btmp.lssolve (atmp, info, rank); |
|
|
227 |
|
1800
|
228 return result.hermitian (); |
|
1
|
229 } |
|
|
230 |
|
767
|
231 // Funny element by element division operations. |
|
|
232 // |
|
|
233 // op2 \ op1: s cs |
|
|
234 // +-- +---+----+ |
|
|
235 // matrix | 1 | 3 | |
|
|
236 // +---+----+ |
|
|
237 // complex_matrix | 2 | 4 | |
|
|
238 // +---+----+ |
|
1
|
239 |
|
1800
|
240 Matrix |
|
164
|
241 x_el_div (double a, const Matrix& b) |
|
1
|
242 { |
|
5275
|
243 octave_idx_type nr = b.rows (); |
|
|
244 octave_idx_type nc = b.columns (); |
|
1
|
245 |
|
|
246 Matrix result (nr, nc); |
|
|
247 |
|
5275
|
248 for (octave_idx_type j = 0; j < nc; j++) |
|
|
249 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
250 { |
|
|
251 OCTAVE_QUIT; |
|
|
252 result (i, j) = a / b (i, j); |
|
|
253 } |
|
1
|
254 |
|
1800
|
255 return result; |
|
1
|
256 } |
|
|
257 |
|
1800
|
258 ComplexMatrix |
|
164
|
259 x_el_div (double a, const ComplexMatrix& b) |
|
1
|
260 { |
|
5275
|
261 octave_idx_type nr = b.rows (); |
|
|
262 octave_idx_type nc = b.columns (); |
|
1
|
263 |
|
|
264 ComplexMatrix result (nr, nc); |
|
|
265 |
|
5275
|
266 for (octave_idx_type j = 0; j < nc; j++) |
|
|
267 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
268 { |
|
|
269 OCTAVE_QUIT; |
|
|
270 result (i, j) = a / b (i, j); |
|
|
271 } |
|
1
|
272 |
|
1800
|
273 return result; |
|
1
|
274 } |
|
|
275 |
|
1800
|
276 ComplexMatrix |
|
164
|
277 x_el_div (const Complex a, const Matrix& b) |
|
1
|
278 { |
|
5275
|
279 octave_idx_type nr = b.rows (); |
|
|
280 octave_idx_type nc = b.columns (); |
|
1
|
281 |
|
|
282 ComplexMatrix result (nr, nc); |
|
|
283 |
|
5275
|
284 for (octave_idx_type j = 0; j < nc; j++) |
|
|
285 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
286 { |
|
|
287 OCTAVE_QUIT; |
|
|
288 result (i, j) = a / b (i, j); |
|
|
289 } |
|
1
|
290 |
|
1800
|
291 return result; |
|
1
|
292 } |
|
|
293 |
|
1800
|
294 ComplexMatrix |
|
164
|
295 x_el_div (const Complex a, const ComplexMatrix& b) |
|
1
|
296 { |
|
5275
|
297 octave_idx_type nr = b.rows (); |
|
|
298 octave_idx_type nc = b.columns (); |
|
1
|
299 |
|
|
300 ComplexMatrix result (nr, nc); |
|
|
301 |
|
5275
|
302 for (octave_idx_type j = 0; j < nc; j++) |
|
|
303 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
304 { |
|
|
305 OCTAVE_QUIT; |
|
|
306 result (i, j) = a / b (i, j); |
|
|
307 } |
|
1
|
308 |
|
1800
|
309 return result; |
|
1
|
310 } |
|
|
311 |
|
4543
|
312 // Funny element by element division operations. |
|
|
313 // |
|
|
314 // op2 \ op1: s cs |
|
|
315 // +-- +---+----+ |
|
|
316 // N-d array | 1 | 3 | |
|
|
317 // +---+----+ |
|
|
318 // complex N-d array | 2 | 4 | |
|
|
319 // +---+----+ |
|
|
320 |
|
|
321 NDArray |
|
|
322 x_el_div (double a, const NDArray& b) |
|
|
323 { |
|
|
324 NDArray result (b.dims ()); |
|
|
325 |
|
5275
|
326 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
327 { |
|
|
328 OCTAVE_QUIT; |
|
|
329 result (i) = a / b (i); |
|
|
330 } |
|
|
331 |
|
|
332 return result; |
|
|
333 } |
|
|
334 |
|
|
335 ComplexNDArray |
|
|
336 x_el_div (double a, const ComplexNDArray& b) |
|
|
337 { |
|
|
338 ComplexNDArray result (b.dims ()); |
|
|
339 |
|
5275
|
340 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
341 { |
|
|
342 OCTAVE_QUIT; |
|
|
343 result (i) = a / b (i); |
|
|
344 } |
|
|
345 |
|
|
346 return result; |
|
|
347 } |
|
|
348 |
|
|
349 ComplexNDArray |
|
|
350 x_el_div (const Complex a, const NDArray& b) |
|
|
351 { |
|
|
352 ComplexNDArray result (b.dims ()); |
|
|
353 |
|
5275
|
354 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
355 { |
|
|
356 OCTAVE_QUIT; |
|
|
357 result (i) = a / b (i); |
|
|
358 } |
|
|
359 |
|
|
360 return result; |
|
|
361 } |
|
|
362 |
|
|
363 ComplexNDArray |
|
|
364 x_el_div (const Complex a, const ComplexNDArray& b) |
|
|
365 { |
|
|
366 ComplexNDArray result (b.dims ()); |
|
|
367 |
|
5275
|
368 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
369 { |
|
|
370 OCTAVE_QUIT; |
|
|
371 result (i) = a / b (i); |
|
|
372 } |
|
|
373 |
|
|
374 return result; |
|
|
375 } |
|
|
376 |
|
767
|
377 // Left division functions. |
|
|
378 // |
|
|
379 // op2 \ op1: m cm |
|
|
380 // +-- +---+----+ |
|
|
381 // matrix | 1 | 3 | |
|
|
382 // +---+----+ |
|
|
383 // complex_matrix | 2 | 4 | |
|
|
384 // +---+----+ |
|
1
|
385 |
|
767
|
386 // -*- 1 -*- |
|
1800
|
387 Matrix |
|
164
|
388 xleftdiv (const Matrix& a, const Matrix& b) |
|
1
|
389 { |
|
2364
|
390 if (! mx_leftdiv_conform (a, b)) |
|
1800
|
391 return Matrix (); |
|
1
|
392 |
|
5275
|
393 octave_idx_type info; |
|
1
|
394 if (a.rows () == a.columns ()) |
|
|
395 { |
|
|
396 double rcond = 0.0; |
|
3480
|
397 |
|
|
398 Matrix result |
|
|
399 = a.solve (b, info, rcond, solve_singularity_warning); |
|
|
400 |
|
|
401 if (result_ok (info)) |
|
1800
|
402 return result; |
|
1
|
403 } |
|
|
404 |
|
5275
|
405 octave_idx_type rank; |
|
1800
|
406 return a.lssolve (b, info, rank); |
|
1
|
407 } |
|
|
408 |
|
767
|
409 // -*- 2 -*- |
|
1800
|
410 ComplexMatrix |
|
164
|
411 xleftdiv (const Matrix& a, const ComplexMatrix& b) |
|
1
|
412 { |
|
2364
|
413 if (! mx_leftdiv_conform (a, b)) |
|
1800
|
414 return ComplexMatrix (); |
|
1
|
415 |
|
5275
|
416 octave_idx_type info; |
|
1
|
417 if (a.rows () == a.columns ()) |
|
|
418 { |
|
|
419 double rcond = 0.0; |
|
3480
|
420 |
|
|
421 ComplexMatrix result |
|
|
422 = a.solve (b, info, rcond, solve_singularity_warning); |
|
|
423 |
|
|
424 if (result_ok (info)) |
|
1800
|
425 return result; |
|
1
|
426 } |
|
|
427 |
|
5275
|
428 octave_idx_type rank; |
|
1800
|
429 return a.lssolve (b, info, rank); |
|
1
|
430 } |
|
|
431 |
|
767
|
432 // -*- 3 -*- |
|
1800
|
433 ComplexMatrix |
|
164
|
434 xleftdiv (const ComplexMatrix& a, const Matrix& b) |
|
1
|
435 { |
|
2364
|
436 if (! mx_leftdiv_conform (a, b)) |
|
1800
|
437 return ComplexMatrix (); |
|
1
|
438 |
|
5275
|
439 octave_idx_type info; |
|
1
|
440 if (a.rows () == a.columns ()) |
|
|
441 { |
|
|
442 double rcond = 0.0; |
|
3480
|
443 |
|
|
444 ComplexMatrix result |
|
|
445 = a.solve (b, info, rcond, solve_singularity_warning); |
|
|
446 |
|
|
447 if (result_ok (info)) |
|
1800
|
448 return result; |
|
1
|
449 } |
|
|
450 |
|
5275
|
451 octave_idx_type rank; |
|
1800
|
452 return a.lssolve (b, info, rank); |
|
1
|
453 } |
|
|
454 |
|
767
|
455 // -*- 4 -*- |
|
1800
|
456 ComplexMatrix |
|
164
|
457 xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b) |
|
1
|
458 { |
|
2364
|
459 if (! mx_leftdiv_conform (a, b)) |
|
1800
|
460 return ComplexMatrix (); |
|
1
|
461 |
|
5275
|
462 octave_idx_type info; |
|
1
|
463 if (a.rows () == a.columns ()) |
|
|
464 { |
|
|
465 double rcond = 0.0; |
|
3480
|
466 |
|
|
467 ComplexMatrix result |
|
|
468 = a.solve (b, info, rcond, solve_singularity_warning); |
|
|
469 |
|
|
470 if (result_ok (info)) |
|
1800
|
471 return result; |
|
1
|
472 } |
|
|
473 |
|
5275
|
474 octave_idx_type rank; |
|
1800
|
475 return a.lssolve (b, info, rank); |
|
1
|
476 } |
|
|
477 |
|
|
478 /* |
|
|
479 ;;; Local Variables: *** |
|
|
480 ;;; mode: C++ *** |
|
|
481 ;;; End: *** |
|
|
482 */ |
