5164
|
1 /* |
|
2 |
|
3 Copyright (C) 2004 David Bateman |
|
4 Copyright (C) 1998-2004 Andy Adler |
|
5 |
|
6 Octave is free software; you can redistribute it and/or modify it |
|
7 under the terms of the GNU General Public License as published by the |
|
8 Free Software Foundation; either version 2, or (at your option) any |
|
9 later version. |
|
10 |
|
11 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
14 for more details. |
|
15 |
|
16 You should have received a copy of the GNU General Public License |
|
17 along with this program; see the file COPYING. If not, write to the Free |
|
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
19 |
|
20 */ |
|
21 |
|
22 #ifdef HAVE_CONFIG_H |
|
23 #include <config.h> |
|
24 #endif |
|
25 |
|
26 #include <cfloat> |
|
27 |
|
28 #include <iostream> |
|
29 #include <vector> |
|
30 |
|
31 #include "quit.h" |
|
32 #include "lo-ieee.h" |
|
33 #include "lo-mappers.h" |
|
34 #include "f77-fcn.h" |
|
35 #include "dRowVector.h" |
|
36 |
|
37 #include "CSparse.h" |
|
38 #include "boolSparse.h" |
|
39 #include "dSparse.h" |
|
40 #include "oct-spparms.h" |
|
41 #include "SparseCmplxLU.h" |
|
42 |
|
43 // External UMFPACK functions in C |
|
44 extern "C" { |
|
45 #include "umfpack.h" |
|
46 } |
|
47 |
|
48 // Fortran functions we call. |
|
49 extern "C" |
|
50 { |
|
51 F77_RET_T |
|
52 F77_FUNC (zgbtrf, ZGBTRF) (const int&, const int&, const int&, |
|
53 const int&, Complex*, const int&, int*, int&); |
|
54 |
|
55 F77_RET_T |
|
56 F77_FUNC (zgbtrs, ZGBTRS) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
57 const int&, const int&, const int&, |
|
58 const Complex*, const int&, |
|
59 const int*, Complex*, const int&, int& |
|
60 F77_CHAR_ARG_LEN_DECL); |
|
61 |
|
62 F77_RET_T |
|
63 F77_FUNC (zgbcon, ZGBCON) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
64 const int&, const int&, Complex*, |
|
65 const int&, const int*, const double&, |
|
66 double&, Complex*, double*, int& |
|
67 F77_CHAR_ARG_LEN_DECL); |
|
68 |
|
69 F77_RET_T |
|
70 F77_FUNC (zpbtrf, ZPBTRF) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
71 const int&, Complex*, const int&, int& |
|
72 F77_CHAR_ARG_LEN_DECL); |
|
73 |
|
74 F77_RET_T |
|
75 F77_FUNC (zpbtrs, ZPBTRS) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
76 const int&, const int&, Complex*, const int&, |
|
77 Complex*, const int&, int& |
|
78 F77_CHAR_ARG_LEN_DECL); |
|
79 |
|
80 F77_RET_T |
|
81 F77_FUNC (zpbcon, ZPBCON) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
82 const int&, Complex*, const int&, |
|
83 const double&, double&, Complex*, int*, int& |
|
84 F77_CHAR_ARG_LEN_DECL); |
|
85 |
|
86 F77_RET_T |
|
87 F77_FUNC (zgttrf, ZGTTRF) (const int&, Complex*, Complex*, Complex*, |
|
88 Complex*, int*, int&); |
|
89 |
|
90 F77_RET_T |
|
91 F77_FUNC (zgttrs, ZGTTRS) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
92 const int&, const Complex*, const Complex*, |
|
93 const Complex*, const Complex*, const int*, |
|
94 Complex *, const int&, int& |
|
95 F77_CHAR_ARG_LEN_DECL); |
|
96 |
|
97 F77_RET_T |
|
98 F77_FUNC (zptsv, ZPTSV) (const int&, const int&, Complex*, Complex*, |
|
99 Complex*, const int&, int&); |
|
100 |
|
101 F77_RET_T |
|
102 F77_FUNC (zgtsv, ZGTSV) (const int&, const int&, Complex*, Complex*, |
|
103 Complex*, Complex*, const int&, int&); |
|
104 } |
|
105 |
|
106 SparseComplexMatrix::SparseComplexMatrix (const SparseMatrix& a) |
|
107 : MSparse<Complex> (a.rows (), a.cols (), a.nnz ()) |
|
108 { |
|
109 int nc = cols (); |
|
110 int nz = nnz (); |
|
111 |
|
112 for (int i = 0; i < nc + 1; i++) |
|
113 cidx (i) = a.cidx (i); |
|
114 |
|
115 for (int i = 0; i < nz; i++) |
|
116 { |
|
117 data (i) = a.data (i); |
|
118 ridx (i) = a.ridx (i); |
|
119 } |
|
120 } |
|
121 |
|
122 SparseComplexMatrix::SparseComplexMatrix (const SparseBoolMatrix& a) |
|
123 : MSparse<Complex> (a.rows (), a.cols (), a.nnz ()) |
|
124 { |
|
125 int nc = cols (); |
|
126 int nz = nnz (); |
|
127 |
|
128 for (int i = 0; i < nc + 1; i++) |
|
129 cidx (i) = a.cidx (i); |
|
130 |
|
131 for (int i = 0; i < nz; i++) |
|
132 { |
|
133 data (i) = a.data (i); |
|
134 ridx (i) = a.ridx (i); |
|
135 } |
|
136 } |
|
137 |
|
138 bool |
|
139 SparseComplexMatrix::operator == (const SparseComplexMatrix& a) const |
|
140 { |
|
141 int nr = rows (); |
|
142 int nc = cols (); |
|
143 int nz = nnz (); |
|
144 int nr_a = a.rows (); |
|
145 int nc_a = a.cols (); |
|
146 int nz_a = a.nnz (); |
|
147 |
|
148 if (nr != nr_a || nc != nc_a || nz != nz_a) |
|
149 return false; |
|
150 |
|
151 for (int i = 0; i < nc + 1; i++) |
|
152 if (cidx(i) != a.cidx(i)) |
|
153 return false; |
|
154 |
|
155 for (int i = 0; i < nz; i++) |
|
156 if (data(i) != a.data(i) || ridx(i) != a.ridx(i)) |
|
157 return false; |
|
158 |
|
159 return true; |
|
160 } |
|
161 |
|
162 bool |
|
163 SparseComplexMatrix::operator != (const SparseComplexMatrix& a) const |
|
164 { |
|
165 return !(*this == a); |
|
166 } |
|
167 |
|
168 bool |
|
169 SparseComplexMatrix::is_hermitian (void) const |
|
170 { |
|
171 int nr = rows (); |
|
172 int nc = cols (); |
|
173 |
|
174 if (is_square () && nr > 0) |
|
175 { |
|
176 for (int i = 0; i < nr; i++) |
|
177 for (int j = i; j < nc; j++) |
|
178 if (elem (i, j) != conj (elem (j, i))) |
|
179 return false; |
|
180 |
|
181 return true; |
|
182 } |
|
183 |
|
184 return false; |
|
185 } |
|
186 |
|
187 static const Complex Complex_NaN_result (octave_NaN, octave_NaN); |
|
188 |
|
189 SparseComplexMatrix |
|
190 SparseComplexMatrix::max (int dim) const |
|
191 { |
|
192 Array2<int> dummy_idx; |
|
193 return max (dummy_idx, dim); |
|
194 } |
|
195 |
|
196 SparseComplexMatrix |
|
197 SparseComplexMatrix::max (Array2<int>& idx_arg, int dim) const |
|
198 { |
|
199 SparseComplexMatrix result; |
|
200 dim_vector dv = dims (); |
|
201 |
|
202 if (dv.numel () == 0 || dim > dv.length () || dim < 0) |
|
203 return result; |
|
204 |
|
205 int nr = dv(0); |
|
206 int nc = dv(1); |
|
207 |
|
208 if (dim == 0) |
|
209 { |
|
210 idx_arg.resize (1, nc); |
|
211 int nel = 0; |
|
212 for (int j = 0; j < nc; j++) |
|
213 { |
|
214 Complex tmp_max; |
|
215 double abs_max = octave_NaN; |
|
216 int idx_j = 0; |
|
217 for (int i = cidx(j); i < cidx(j+1); i++) |
|
218 { |
|
219 if (ridx(i) != idx_j) |
|
220 break; |
|
221 else |
|
222 idx_j++; |
|
223 } |
|
224 |
|
225 if (idx_j != nr) |
|
226 { |
|
227 tmp_max = 0.; |
|
228 abs_max = 0.; |
|
229 } |
|
230 |
|
231 for (int i = cidx(j); i < cidx(j+1); i++) |
|
232 { |
|
233 Complex tmp = data (i); |
|
234 |
|
235 if (octave_is_NaN_or_NA (tmp)) |
|
236 continue; |
|
237 |
|
238 double abs_tmp = ::abs (tmp); |
|
239 |
|
240 if (octave_is_NaN_or_NA (abs_max) || abs_tmp > abs_max) |
|
241 { |
|
242 idx_j = ridx (i); |
|
243 tmp_max = tmp; |
|
244 abs_max = abs_tmp; |
|
245 } |
|
246 } |
|
247 |
|
248 idx_arg.elem (j) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_j; |
|
249 if (abs_max != 0.) |
|
250 nel++; |
|
251 } |
|
252 |
|
253 result = SparseComplexMatrix (1, nc, nel); |
|
254 |
|
255 int ii = 0; |
|
256 result.xcidx (0) = 0; |
|
257 for (int j = 0; j < nc; j++) |
|
258 { |
|
259 Complex tmp = elem (idx_arg(j), j); |
|
260 if (tmp != 0.) |
|
261 { |
|
262 result.xdata (ii) = tmp; |
|
263 result.xridx (ii++) = 0; |
|
264 } |
|
265 result.xcidx (j+1) = ii; |
|
266 } |
|
267 } |
|
268 else |
|
269 { |
|
270 idx_arg.resize (nr, 1, 0); |
|
271 |
|
272 for (int i = cidx(0); i < cidx(1); i++) |
|
273 idx_arg.elem(ridx(i)) = -1; |
|
274 |
|
275 for (int j = 0; j < nc; j++) |
|
276 for (int i = 0; i < nr; i++) |
|
277 { |
|
278 if (idx_arg.elem(i) != -1) |
|
279 continue; |
|
280 bool found = false; |
|
281 for (int k = cidx(j); k < cidx(j+1); k++) |
|
282 if (ridx(k) == i) |
|
283 { |
|
284 found = true; |
|
285 break; |
|
286 } |
|
287 |
|
288 if (!found) |
|
289 idx_arg.elem(i) = j; |
|
290 |
|
291 } |
|
292 |
|
293 for (int j = 0; j < nc; j++) |
|
294 { |
|
295 for (int i = cidx(j); i < cidx(j+1); i++) |
|
296 { |
|
297 int ir = ridx (i); |
|
298 int ix = idx_arg.elem (ir); |
|
299 Complex tmp = data (i); |
|
300 |
|
301 if (octave_is_NaN_or_NA (tmp)) |
|
302 continue; |
|
303 else if (ix == -1 || ::abs(tmp) > ::abs(elem (ir, ix))) |
|
304 idx_arg.elem (ir) = j; |
|
305 } |
|
306 } |
|
307 |
|
308 int nel = 0; |
|
309 for (int j = 0; j < nr; j++) |
|
310 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) |
|
311 nel++; |
|
312 |
|
313 result = SparseComplexMatrix (nr, 1, nel); |
|
314 |
|
315 int ii = 0; |
|
316 result.xcidx (0) = 0; |
|
317 result.xcidx (1) = nel; |
|
318 for (int j = 0; j < nr; j++) |
|
319 { |
|
320 if (idx_arg(j) == -1) |
|
321 { |
|
322 idx_arg(j) = 0; |
|
323 result.xdata (ii) = Complex_NaN_result; |
|
324 result.xridx (ii++) = j; |
|
325 } |
|
326 else |
|
327 { |
|
328 Complex tmp = elem (j, idx_arg(j)); |
|
329 if (tmp != 0.) |
|
330 { |
|
331 result.xdata (ii) = tmp; |
|
332 result.xridx (ii++) = j; |
|
333 } |
|
334 } |
|
335 } |
|
336 } |
|
337 |
|
338 return result; |
|
339 } |
|
340 |
|
341 SparseComplexMatrix |
|
342 SparseComplexMatrix::min (int dim) const |
|
343 { |
|
344 Array2<int> dummy_idx; |
|
345 return min (dummy_idx, dim); |
|
346 } |
|
347 |
|
348 SparseComplexMatrix |
|
349 SparseComplexMatrix::min (Array2<int>& idx_arg, int dim) const |
|
350 { |
|
351 SparseComplexMatrix result; |
|
352 dim_vector dv = dims (); |
|
353 |
|
354 if (dv.numel () == 0 || dim > dv.length () || dim < 0) |
|
355 return result; |
|
356 |
|
357 int nr = dv(0); |
|
358 int nc = dv(1); |
|
359 |
|
360 if (dim == 0) |
|
361 { |
|
362 idx_arg.resize (1, nc); |
|
363 int nel = 0; |
|
364 for (int j = 0; j < nc; j++) |
|
365 { |
|
366 Complex tmp_min; |
|
367 double abs_min = octave_NaN; |
|
368 int idx_j = 0; |
|
369 for (int i = cidx(j); i < cidx(j+1); i++) |
|
370 { |
|
371 if (ridx(i) != idx_j) |
|
372 break; |
|
373 else |
|
374 idx_j++; |
|
375 } |
|
376 |
|
377 if (idx_j != nr) |
|
378 { |
|
379 tmp_min = 0.; |
|
380 abs_min = 0.; |
|
381 } |
|
382 |
|
383 for (int i = cidx(j); i < cidx(j+1); i++) |
|
384 { |
|
385 Complex tmp = data (i); |
|
386 |
|
387 if (octave_is_NaN_or_NA (tmp)) |
|
388 continue; |
|
389 |
|
390 double abs_tmp = ::abs (tmp); |
|
391 |
|
392 if (octave_is_NaN_or_NA (abs_min) || abs_tmp < abs_min) |
|
393 { |
|
394 idx_j = ridx (i); |
|
395 tmp_min = tmp; |
|
396 abs_min = abs_tmp; |
|
397 } |
|
398 } |
|
399 |
|
400 idx_arg.elem (j) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_j; |
|
401 if (abs_min != 0.) |
|
402 nel++; |
|
403 } |
|
404 |
|
405 result = SparseComplexMatrix (1, nc, nel); |
|
406 |
|
407 int ii = 0; |
|
408 result.xcidx (0) = 0; |
|
409 for (int j = 0; j < nc; j++) |
|
410 { |
|
411 Complex tmp = elem (idx_arg(j), j); |
|
412 if (tmp != 0.) |
|
413 { |
|
414 result.xdata (ii) = tmp; |
|
415 result.xridx (ii++) = 0; |
|
416 } |
|
417 result.xcidx (j+1) = ii; |
|
418 } |
|
419 } |
|
420 else |
|
421 { |
|
422 idx_arg.resize (nr, 1, 0); |
|
423 |
|
424 for (int i = cidx(0); i < cidx(1); i++) |
|
425 idx_arg.elem(ridx(i)) = -1; |
|
426 |
|
427 for (int j = 0; j < nc; j++) |
|
428 for (int i = 0; i < nr; i++) |
|
429 { |
|
430 if (idx_arg.elem(i) != -1) |
|
431 continue; |
|
432 bool found = false; |
|
433 for (int k = cidx(j); k < cidx(j+1); k++) |
|
434 if (ridx(k) == i) |
|
435 { |
|
436 found = true; |
|
437 break; |
|
438 } |
|
439 |
|
440 if (!found) |
|
441 idx_arg.elem(i) = j; |
|
442 |
|
443 } |
|
444 |
|
445 for (int j = 0; j < nc; j++) |
|
446 { |
|
447 for (int i = cidx(j); i < cidx(j+1); i++) |
|
448 { |
|
449 int ir = ridx (i); |
|
450 int ix = idx_arg.elem (ir); |
|
451 Complex tmp = data (i); |
|
452 |
|
453 if (octave_is_NaN_or_NA (tmp)) |
|
454 continue; |
|
455 else if (ix == -1 || ::abs(tmp) < ::abs(elem (ir, ix))) |
|
456 idx_arg.elem (ir) = j; |
|
457 } |
|
458 } |
|
459 |
|
460 int nel = 0; |
|
461 for (int j = 0; j < nr; j++) |
|
462 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) |
|
463 nel++; |
|
464 |
|
465 result = SparseComplexMatrix (nr, 1, nel); |
|
466 |
|
467 int ii = 0; |
|
468 result.xcidx (0) = 0; |
|
469 result.xcidx (1) = nel; |
|
470 for (int j = 0; j < nr; j++) |
|
471 { |
|
472 if (idx_arg(j) == -1) |
|
473 { |
|
474 idx_arg(j) = 0; |
|
475 result.xdata (ii) = Complex_NaN_result; |
|
476 result.xridx (ii++) = j; |
|
477 } |
|
478 else |
|
479 { |
|
480 Complex tmp = elem (j, idx_arg(j)); |
|
481 if (tmp != 0.) |
|
482 { |
|
483 result.xdata (ii) = tmp; |
|
484 result.xridx (ii++) = j; |
|
485 } |
|
486 } |
|
487 } |
|
488 } |
|
489 |
|
490 return result; |
|
491 } |
|
492 |
|
493 // destructive insert/delete/reorder operations |
|
494 |
|
495 SparseComplexMatrix& |
|
496 SparseComplexMatrix::insert (const SparseMatrix& a, int r, int c) |
|
497 { |
|
498 SparseComplexMatrix tmp (a); |
|
499 return insert (a, r, c); |
|
500 } |
|
501 |
|
502 SparseComplexMatrix& |
|
503 SparseComplexMatrix::insert (const SparseComplexMatrix& a, int r, int c) |
|
504 { |
|
505 MSparse<Complex>::insert (a, r, c); |
|
506 return *this; |
|
507 } |
|
508 |
|
509 SparseComplexMatrix |
|
510 SparseComplexMatrix::concat (const SparseComplexMatrix& rb, |
|
511 const Array<int>& ra_idx) |
|
512 { |
|
513 // Don't use numel to avoid all possiblity of an overflow |
|
514 if (rb.rows () > 0 && rb.cols () > 0) |
|
515 insert (rb, ra_idx(0), ra_idx(1)); |
|
516 return *this; |
|
517 } |
|
518 |
|
519 SparseComplexMatrix |
|
520 SparseComplexMatrix::concat (const SparseMatrix& rb, const Array<int>& ra_idx) |
|
521 { |
|
522 SparseComplexMatrix tmp (rb); |
|
523 if (rb.rows () > 0 && rb.cols () > 0) |
|
524 insert (tmp, ra_idx(0), ra_idx(1)); |
|
525 return *this; |
|
526 } |
|
527 |
|
528 ComplexMatrix |
|
529 SparseComplexMatrix::matrix_value (void) const |
|
530 { |
|
531 int nr = rows (); |
|
532 int nc = cols (); |
|
533 ComplexMatrix retval (nr, nc, Complex (0.0, 0.0)); |
|
534 |
|
535 for (int j = 0; j < nc; j++) |
|
536 for (int i = cidx(j); i < cidx(j+1); i++) |
|
537 retval.elem (ridx(i), j) = data (i); |
|
538 |
|
539 return retval; |
|
540 } |
|
541 |
|
542 SparseComplexMatrix |
|
543 SparseComplexMatrix::hermitian (void) const |
|
544 { |
|
545 int nr = rows (); |
|
546 int nc = cols (); |
|
547 int nz = nnz (); |
|
548 SparseComplexMatrix retval (nc, nr, nz); |
|
549 |
|
550 retval.cidx(0) = 0; |
|
551 for (int i = 0, iidx = 0; i < nr; i++) |
|
552 { |
|
553 for (int j = 0; j < nc; j++) |
|
554 for (int k = cidx(j); k < cidx(j+1); k++) |
|
555 if (ridx(k) == i) |
|
556 { |
|
557 retval.data(iidx) = conj (data(k)); |
|
558 retval.ridx(iidx++) = j; |
|
559 } |
|
560 retval.cidx(i+1) = iidx; |
|
561 } |
|
562 |
|
563 return retval; |
|
564 } |
|
565 |
|
566 SparseComplexMatrix |
|
567 conj (const SparseComplexMatrix& a) |
|
568 { |
|
569 int nr = a.rows (); |
|
570 int nc = a.cols (); |
|
571 int nz = a.nnz (); |
|
572 SparseComplexMatrix retval (nc, nr, nz); |
|
573 |
|
574 for (int i = 0; i < nc + 1; i++) |
|
575 retval.cidx (i) = a.cidx (i); |
|
576 |
|
577 for (int i = 0; i < nz; i++) |
|
578 { |
|
579 retval.data (i) = conj (a.data (i)); |
|
580 retval.ridx (i) = a.ridx (i); |
|
581 } |
|
582 |
|
583 return retval; |
|
584 } |
|
585 |
|
586 SparseComplexMatrix |
|
587 SparseComplexMatrix::inverse (void) const |
|
588 { |
|
589 int info; |
|
590 double rcond; |
|
591 return inverse (info, rcond, 0, 0); |
|
592 } |
|
593 |
|
594 SparseComplexMatrix |
|
595 SparseComplexMatrix::inverse (int& info) const |
|
596 { |
|
597 double rcond; |
|
598 return inverse (info, rcond, 0, 0); |
|
599 } |
|
600 |
|
601 SparseComplexMatrix |
|
602 SparseComplexMatrix::inverse (int& info, double& rcond, int force, |
|
603 int calc_cond) const |
|
604 { |
|
605 info = -1; |
|
606 (*current_liboctave_error_handler) |
|
607 ("SparseComplexMatrix::inverse not implemented yet"); |
|
608 return SparseComplexMatrix (); |
|
609 } |
|
610 |
|
611 ComplexDET |
|
612 SparseComplexMatrix::determinant (void) const |
|
613 { |
|
614 int info; |
|
615 double rcond; |
|
616 return determinant (info, rcond, 0); |
|
617 } |
|
618 |
|
619 ComplexDET |
|
620 SparseComplexMatrix::determinant (int& info) const |
|
621 { |
|
622 double rcond; |
|
623 return determinant (info, rcond, 0); |
|
624 } |
|
625 |
|
626 ComplexDET |
|
627 SparseComplexMatrix::determinant (int& err, double& rcond, int calc_cond) const |
|
628 { |
|
629 ComplexDET retval; |
|
630 |
|
631 int nr = rows (); |
|
632 int nc = cols (); |
|
633 |
|
634 if (nr == 0 || nc == 0 || nr != nc) |
|
635 { |
|
636 Complex d[2]; |
|
637 d[0] = 1.0; |
|
638 d[1] = 0.0; |
|
639 retval = ComplexDET (d); |
|
640 } |
|
641 else |
|
642 { |
|
643 err = 0; |
|
644 |
|
645 // Setup the control parameters |
|
646 Matrix Control (UMFPACK_CONTROL, 1); |
|
647 double *control = Control.fortran_vec (); |
|
648 umfpack_zi_defaults (control); |
|
649 |
|
650 double tmp = Voctave_sparse_controls.get_key ("spumoni"); |
|
651 if (!xisnan (tmp)) |
|
652 Control (UMFPACK_PRL) = tmp; |
|
653 |
|
654 tmp = Voctave_sparse_controls.get_key ("piv_tol"); |
|
655 if (!xisnan (tmp)) |
|
656 { |
|
657 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
|
658 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
659 } |
|
660 |
|
661 // Set whether we are allowed to modify Q or not |
|
662 tmp = Voctave_sparse_controls.get_key ("autoamd"); |
|
663 if (!xisnan (tmp)) |
|
664 Control (UMFPACK_FIXQ) = tmp; |
|
665 |
|
666 // Turn-off UMFPACK scaling for LU |
|
667 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; |
|
668 |
|
669 umfpack_zi_report_control (control); |
|
670 |
|
671 const int *Ap = cidx (); |
|
672 const int *Ai = ridx (); |
|
673 const Complex *Ax = data (); |
|
674 |
|
675 umfpack_zi_report_matrix (nr, nc, Ap, Ai, |
|
676 X_CAST (const double *, Ax), |
|
677 NULL, 1, control); |
|
678 |
|
679 void *Symbolic; |
|
680 Matrix Info (1, UMFPACK_INFO); |
|
681 double *info = Info.fortran_vec (); |
|
682 int status = umfpack_zi_qsymbolic |
|
683 (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, |
|
684 NULL, &Symbolic, control, info); |
|
685 |
|
686 if (status < 0) |
|
687 { |
|
688 (*current_liboctave_error_handler) |
|
689 ("SparseComplexMatrix::determinant symbolic factorization failed"); |
|
690 |
|
691 umfpack_zi_report_status (control, status); |
|
692 umfpack_zi_report_info (control, info); |
|
693 |
|
694 umfpack_zi_free_symbolic (&Symbolic) ; |
|
695 } |
|
696 else |
|
697 { |
|
698 umfpack_zi_report_symbolic (Symbolic, control); |
|
699 |
|
700 void *Numeric; |
|
701 status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), |
|
702 NULL, Symbolic, &Numeric, |
|
703 control, info) ; |
|
704 umfpack_zi_free_symbolic (&Symbolic) ; |
|
705 |
|
706 rcond = Info (UMFPACK_RCOND); |
|
707 |
|
708 if (status < 0) |
|
709 { |
|
710 (*current_liboctave_error_handler) |
|
711 ("SparseComplexMatrix::determinant numeric factorization failed"); |
|
712 |
|
713 umfpack_zi_report_status (control, status); |
|
714 umfpack_zi_report_info (control, info); |
|
715 |
|
716 umfpack_zi_free_numeric (&Numeric); |
|
717 } |
|
718 else |
|
719 { |
|
720 umfpack_zi_report_numeric (Numeric, control); |
|
721 |
|
722 Complex d[2]; |
|
723 double d_exponent; |
|
724 |
|
725 status = umfpack_zi_get_determinant |
|
726 (X_CAST (double *, &d[0]), NULL, &d_exponent, |
|
727 Numeric, info); |
|
728 d[1] = d_exponent; |
|
729 |
|
730 if (status < 0) |
|
731 { |
|
732 (*current_liboctave_error_handler) |
|
733 ("SparseComplexMatrix::determinant error calculating determinant"); |
|
734 |
|
735 umfpack_zi_report_status (control, status); |
|
736 umfpack_zi_report_info (control, info); |
|
737 |
|
738 umfpack_zi_free_numeric (&Numeric); |
|
739 } |
|
740 else |
|
741 retval = ComplexDET (d); |
|
742 } |
|
743 } |
|
744 } |
|
745 |
|
746 return retval; |
|
747 } |
|
748 |
|
749 ComplexMatrix |
|
750 SparseComplexMatrix::dsolve (SparseType &mattype, const Matrix& b, int& err, |
|
751 double& rcond, solve_singularity_handler) const |
|
752 { |
|
753 ComplexMatrix retval; |
|
754 |
|
755 int nr = rows (); |
|
756 int nc = cols (); |
|
757 err = 0; |
|
758 |
|
759 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
760 (*current_liboctave_error_handler) |
|
761 ("matrix dimension mismatch solution of linear equations"); |
|
762 else |
|
763 { |
|
764 // Print spparms("spumoni") info if requested |
|
765 int typ = mattype.type (); |
|
766 mattype.info (); |
|
767 |
|
768 if (typ == SparseType::Diagonal || |
|
769 typ == SparseType::Permuted_Diagonal) |
|
770 { |
|
771 retval.resize (b.rows (), b.cols()); |
|
772 if (typ == SparseType::Diagonal) |
|
773 for (int j = 0; j < b.cols(); j++) |
|
774 for (int i = 0; i < nr; i++) |
|
775 retval(i,j) = b(i,j) / data (i); |
|
776 else |
|
777 for (int j = 0; j < b.cols(); j++) |
|
778 for (int i = 0; i < nr; i++) |
|
779 retval(i,j) = b(ridx(i),j) / data (i); |
|
780 |
|
781 double dmax = 0., dmin = octave_Inf; |
|
782 for (int i = 0; i < nr; i++) |
|
783 { |
|
784 double tmp = ::abs(data(i)); |
|
785 if (tmp > dmax) |
|
786 dmax = tmp; |
|
787 if (tmp < dmin) |
|
788 dmin = tmp; |
|
789 } |
|
790 rcond = dmin / dmax; |
|
791 } |
|
792 else |
|
793 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
794 } |
|
795 |
|
796 return retval; |
|
797 } |
|
798 |
|
799 SparseComplexMatrix |
|
800 SparseComplexMatrix::dsolve (SparseType &mattype, const SparseMatrix& b, |
|
801 int& err, double& rcond, solve_singularity_handler) const |
|
802 { |
|
803 SparseComplexMatrix retval; |
|
804 |
|
805 int nr = rows (); |
|
806 int nc = cols (); |
|
807 err = 0; |
|
808 |
|
809 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
810 (*current_liboctave_error_handler) |
|
811 ("matrix dimension mismatch solution of linear equations"); |
|
812 else |
|
813 { |
|
814 // Print spparms("spumoni") info if requested |
|
815 int typ = mattype.type (); |
|
816 mattype.info (); |
|
817 |
|
818 if (typ == SparseType::Diagonal || |
|
819 typ == SparseType::Permuted_Diagonal) |
|
820 { |
|
821 int b_nr = b.rows (); |
|
822 int b_nc = b.cols (); |
|
823 int b_nz = b.nnz (); |
|
824 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
825 |
|
826 retval.xcidx(0) = 0; |
|
827 int ii = 0; |
|
828 if (typ == SparseType::Diagonal) |
|
829 for (int j = 0; j < b.cols(); j++) |
|
830 { |
|
831 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
832 { |
|
833 retval.xridx (ii) = b.ridx(i); |
|
834 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); |
|
835 } |
|
836 retval.xcidx(j+1) = ii; |
|
837 } |
|
838 else |
|
839 for (int j = 0; j < b.cols(); j++) |
|
840 { |
|
841 for (int i = 0; i < nr; i++) |
|
842 { |
|
843 bool found = false; |
|
844 int k; |
|
845 for (k = b.cidx(j); k < b.cidx(j+1); k++) |
|
846 if (ridx(i) == b.ridx(k)) |
|
847 { |
|
848 found = true; |
|
849 break; |
|
850 } |
|
851 if (found) |
|
852 { |
|
853 retval.xridx (ii) = i; |
|
854 retval.xdata (ii++) = b.data(k) / data (i); |
|
855 } |
|
856 } |
|
857 retval.xcidx(j+1) = ii; |
|
858 } |
|
859 |
|
860 double dmax = 0., dmin = octave_Inf; |
|
861 for (int i = 0; i < nr; i++) |
|
862 { |
|
863 double tmp = ::abs(data(i)); |
|
864 if (tmp > dmax) |
|
865 dmax = tmp; |
|
866 if (tmp < dmin) |
|
867 dmin = tmp; |
|
868 } |
|
869 rcond = dmin / dmax; |
|
870 } |
|
871 else |
|
872 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
873 } |
|
874 |
|
875 return retval; |
|
876 } |
|
877 |
|
878 ComplexMatrix |
|
879 SparseComplexMatrix::dsolve (SparseType &mattype, const ComplexMatrix& b, |
|
880 int& err, double& rcond, solve_singularity_handler) const |
|
881 { |
|
882 ComplexMatrix retval; |
|
883 |
|
884 int nr = rows (); |
|
885 int nc = cols (); |
|
886 err = 0; |
|
887 |
|
888 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
889 (*current_liboctave_error_handler) |
|
890 ("matrix dimension mismatch solution of linear equations"); |
|
891 else |
|
892 { |
|
893 // Print spparms("spumoni") info if requested |
|
894 int typ = mattype.type (); |
|
895 mattype.info (); |
|
896 |
|
897 if (typ == SparseType::Diagonal || |
|
898 typ == SparseType::Permuted_Diagonal) |
|
899 { |
|
900 retval.resize (b.rows (), b.cols()); |
|
901 if (typ == SparseType::Diagonal) |
|
902 for (int j = 0; j < b.cols(); j++) |
|
903 for (int i = 0; i < nr; i++) |
|
904 retval(i,j) = b(i,j) / data (i); |
|
905 else |
|
906 for (int j = 0; j < b.cols(); j++) |
|
907 for (int i = 0; i < nr; i++) |
|
908 retval(i,j) = b(ridx(i),j) / data (i); |
|
909 |
|
910 double dmax = 0., dmin = octave_Inf; |
|
911 for (int i = 0; i < nr; i++) |
|
912 { |
|
913 double tmp = ::abs(data(i)); |
|
914 if (tmp > dmax) |
|
915 dmax = tmp; |
|
916 if (tmp < dmin) |
|
917 dmin = tmp; |
|
918 } |
|
919 rcond = dmin / dmax; |
|
920 } |
|
921 else |
|
922 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
923 } |
|
924 |
|
925 return retval; |
|
926 } |
|
927 |
|
928 SparseComplexMatrix |
|
929 SparseComplexMatrix::dsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
930 int& err, double& rcond, |
|
931 solve_singularity_handler) const |
|
932 { |
|
933 SparseComplexMatrix retval; |
|
934 |
|
935 int nr = rows (); |
|
936 int nc = cols (); |
|
937 err = 0; |
|
938 |
|
939 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
940 (*current_liboctave_error_handler) |
|
941 ("matrix dimension mismatch solution of linear equations"); |
|
942 else |
|
943 { |
|
944 // Print spparms("spumoni") info if requested |
|
945 int typ = mattype.type (); |
|
946 mattype.info (); |
|
947 |
|
948 if (typ == SparseType::Diagonal || |
|
949 typ == SparseType::Permuted_Diagonal) |
|
950 { |
|
951 int b_nr = b.rows (); |
|
952 int b_nc = b.cols (); |
|
953 int b_nz = b.nnz (); |
|
954 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
955 |
|
956 retval.xcidx(0) = 0; |
|
957 int ii = 0; |
|
958 if (typ == SparseType::Diagonal) |
|
959 for (int j = 0; j < b.cols(); j++) |
|
960 { |
|
961 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
962 { |
|
963 retval.xridx (ii) = b.ridx(i); |
|
964 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); |
|
965 } |
|
966 retval.xcidx(j+1) = ii; |
|
967 } |
|
968 else |
|
969 for (int j = 0; j < b.cols(); j++) |
|
970 { |
|
971 for (int i = 0; i < nr; i++) |
|
972 { |
|
973 bool found = false; |
|
974 int k; |
|
975 for (k = b.cidx(j); k < b.cidx(j+1); k++) |
|
976 if (ridx(i) == b.ridx(k)) |
|
977 { |
|
978 found = true; |
|
979 break; |
|
980 } |
|
981 if (found) |
|
982 { |
|
983 retval.xridx (ii) = i; |
|
984 retval.xdata (ii++) = b.data(k) / data (i); |
|
985 } |
|
986 } |
|
987 retval.xcidx(j+1) = ii; |
|
988 } |
|
989 |
|
990 double dmax = 0., dmin = octave_Inf; |
|
991 for (int i = 0; i < nr; i++) |
|
992 { |
|
993 double tmp = ::abs(data(i)); |
|
994 if (tmp > dmax) |
|
995 dmax = tmp; |
|
996 if (tmp < dmin) |
|
997 dmin = tmp; |
|
998 } |
|
999 rcond = dmin / dmax; |
|
1000 } |
|
1001 else |
|
1002 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1003 } |
|
1004 |
|
1005 return retval; |
|
1006 } |
|
1007 |
|
1008 ComplexMatrix |
|
1009 SparseComplexMatrix::utsolve (SparseType &mattype, const Matrix& b, int& err, |
|
1010 double& rcond, |
|
1011 solve_singularity_handler sing_handler) const |
|
1012 { |
|
1013 ComplexMatrix retval; |
|
1014 |
|
1015 int nr = rows (); |
|
1016 int nc = cols (); |
|
1017 err = 0; |
|
1018 |
|
1019 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1020 (*current_liboctave_error_handler) |
|
1021 ("matrix dimension mismatch solution of linear equations"); |
|
1022 else |
|
1023 { |
|
1024 // Print spparms("spumoni") info if requested |
|
1025 int typ = mattype.type (); |
|
1026 mattype.info (); |
|
1027 |
|
1028 if (typ == SparseType::Permuted_Upper || |
|
1029 typ == SparseType::Upper) |
|
1030 { |
|
1031 double anorm = 0.; |
|
1032 double ainvnorm = 0.; |
|
1033 int b_cols = b.cols (); |
|
1034 rcond = 0.; |
|
1035 |
|
1036 // Calculate the 1-norm of matrix for rcond calculation |
|
1037 for (int j = 0; j < nr; j++) |
|
1038 { |
|
1039 double atmp = 0.; |
|
1040 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1041 atmp += ::abs(data(i)); |
|
1042 if (atmp > anorm) |
|
1043 anorm = atmp; |
|
1044 } |
|
1045 |
|
1046 if (typ == SparseType::Permuted_Upper) |
|
1047 { |
|
1048 retval.resize (b.rows (), b.cols ()); |
|
1049 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1050 int *p_perm = mattype.triangular_row_perm (); |
|
1051 int *q_perm = mattype.triangular_col_perm (); |
|
1052 |
|
1053 (*current_liboctave_warning_handler) |
|
1054 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1055 |
|
1056 for (int j = 0; j < b_cols; j++) |
|
1057 { |
|
1058 for (int i = 0; i < nr; i++) |
|
1059 work[i] = b(i,j); |
|
1060 |
|
1061 for (int k = nr-1; k >= 0; k--) |
|
1062 { |
|
1063 int iidx = q_perm[k]; |
|
1064 if (work[iidx] != 0.) |
|
1065 { |
|
1066 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1067 { |
|
1068 err = -2; |
|
1069 goto triangular_error; |
|
1070 } |
|
1071 |
|
1072 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1073 work[iidx] = tmp; |
|
1074 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1075 { |
|
1076 int idx2 = q_perm[ridx(i)]; |
|
1077 work[idx2] = |
|
1078 work[idx2] - tmp * data(i); |
|
1079 } |
|
1080 } |
|
1081 } |
|
1082 |
|
1083 for (int i = 0; i < nr; i++) |
|
1084 retval (i, j) = work[p_perm[i]]; |
|
1085 } |
|
1086 |
|
1087 // Calculation of 1-norm of inv(*this) |
|
1088 for (int i = 0; i < nr; i++) |
|
1089 work[i] = 0.; |
|
1090 |
|
1091 for (int j = 0; j < nr; j++) |
|
1092 { |
|
1093 work[q_perm[j]] = 1.; |
|
1094 |
|
1095 for (int k = j; k >= 0; k--) |
|
1096 { |
|
1097 int iidx = q_perm[k]; |
|
1098 |
|
1099 if (work[iidx] != 0.) |
|
1100 { |
|
1101 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1102 work[iidx] = tmp; |
|
1103 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1104 { |
|
1105 int idx2 = q_perm[ridx(i)]; |
|
1106 work[idx2] = work[idx2] - tmp * data(i); |
|
1107 } |
|
1108 } |
|
1109 } |
|
1110 double atmp = 0; |
|
1111 for (int i = 0; i < j+1; i++) |
|
1112 { |
|
1113 atmp += ::abs(work[i]); |
|
1114 work[i] = 0.; |
|
1115 } |
|
1116 if (atmp > ainvnorm) |
|
1117 ainvnorm = atmp; |
|
1118 } |
|
1119 } |
|
1120 else |
|
1121 { |
|
1122 retval = ComplexMatrix (b); |
|
1123 Complex *x_vec = retval.fortran_vec (); |
|
1124 |
|
1125 for (int j = 0; j < b_cols; j++) |
|
1126 { |
|
1127 int offset = j * nr; |
|
1128 for (int k = nr-1; k >= 0; k--) |
|
1129 { |
|
1130 if (x_vec[k+offset] != 0.) |
|
1131 { |
|
1132 if (ridx(cidx(k+1)-1) != k) |
|
1133 { |
|
1134 err = -2; |
|
1135 goto triangular_error; |
|
1136 } |
|
1137 |
|
1138 Complex tmp = x_vec[k+offset] / |
|
1139 data(cidx(k+1)-1); |
|
1140 x_vec[k+offset] = tmp; |
|
1141 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1142 { |
|
1143 int iidx = ridx(i); |
|
1144 x_vec[iidx+offset] = |
|
1145 x_vec[iidx+offset] - tmp * data(i); |
|
1146 } |
|
1147 } |
|
1148 } |
|
1149 } |
|
1150 |
|
1151 // Calculation of 1-norm of inv(*this) |
|
1152 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1153 for (int i = 0; i < nr; i++) |
|
1154 work[i] = 0.; |
|
1155 |
|
1156 for (int j = 0; j < nr; j++) |
|
1157 { |
|
1158 work[j] = 1.; |
|
1159 |
|
1160 for (int k = j; k >= 0; k--) |
|
1161 { |
|
1162 if (work[k] != 0.) |
|
1163 { |
|
1164 Complex tmp = work[k] / data(cidx(k+1)-1); |
|
1165 work[k] = tmp; |
|
1166 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1167 { |
|
1168 int iidx = ridx(i); |
|
1169 work[iidx] = work[iidx] - tmp * data(i); |
|
1170 } |
|
1171 } |
|
1172 } |
|
1173 double atmp = 0; |
|
1174 for (int i = 0; i < j+1; i++) |
|
1175 { |
|
1176 atmp += ::abs(work[i]); |
|
1177 work[i] = 0.; |
|
1178 } |
|
1179 if (atmp > ainvnorm) |
|
1180 ainvnorm = atmp; |
|
1181 } |
|
1182 } |
|
1183 |
|
1184 rcond = 1. / ainvnorm / anorm; |
|
1185 |
|
1186 triangular_error: |
|
1187 if (err != 0) |
|
1188 { |
|
1189 if (sing_handler) |
|
1190 sing_handler (rcond); |
|
1191 else |
|
1192 (*current_liboctave_error_handler) |
|
1193 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
1194 rcond); |
|
1195 } |
|
1196 |
|
1197 volatile double rcond_plus_one = rcond + 1.0; |
|
1198 |
|
1199 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
1200 { |
|
1201 err = -2; |
|
1202 |
|
1203 if (sing_handler) |
|
1204 sing_handler (rcond); |
|
1205 else |
|
1206 (*current_liboctave_error_handler) |
|
1207 ("matrix singular to machine precision, rcond = %g", |
|
1208 rcond); |
|
1209 } |
|
1210 } |
|
1211 else |
|
1212 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1213 } |
|
1214 |
|
1215 return retval; |
|
1216 } |
|
1217 |
|
1218 SparseComplexMatrix |
|
1219 SparseComplexMatrix::utsolve (SparseType &mattype, const SparseMatrix& b, |
|
1220 int& err, double& rcond, |
|
1221 solve_singularity_handler sing_handler) const |
|
1222 { |
|
1223 SparseComplexMatrix retval; |
|
1224 |
|
1225 int nr = rows (); |
|
1226 int nc = cols (); |
|
1227 err = 0; |
|
1228 |
|
1229 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1230 (*current_liboctave_error_handler) |
|
1231 ("matrix dimension mismatch solution of linear equations"); |
|
1232 else |
|
1233 { |
|
1234 // Print spparms("spumoni") info if requested |
|
1235 int typ = mattype.type (); |
|
1236 mattype.info (); |
|
1237 |
|
1238 if (typ == SparseType::Permuted_Upper || |
|
1239 typ == SparseType::Upper) |
|
1240 { |
|
1241 double anorm = 0.; |
|
1242 double ainvnorm = 0.; |
|
1243 rcond = 0.; |
|
1244 |
|
1245 // Calculate the 1-norm of matrix for rcond calculation |
|
1246 for (int j = 0; j < nr; j++) |
|
1247 { |
|
1248 double atmp = 0.; |
|
1249 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1250 atmp += ::abs(data(i)); |
|
1251 if (atmp > anorm) |
|
1252 anorm = atmp; |
|
1253 } |
|
1254 |
|
1255 int b_nr = b.rows (); |
|
1256 int b_nc = b.cols (); |
|
1257 int b_nz = b.nnz (); |
|
1258 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
1259 retval.xcidx(0) = 0; |
|
1260 int ii = 0; |
|
1261 int x_nz = b_nz; |
|
1262 |
|
1263 if (typ == SparseType::Permuted_Upper) |
|
1264 { |
|
1265 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1266 int *p_perm = mattype.triangular_row_perm (); |
|
1267 int *q_perm = mattype.triangular_col_perm (); |
|
1268 |
|
1269 (*current_liboctave_warning_handler) |
|
1270 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1271 |
|
1272 for (int j = 0; j < b_nc; j++) |
|
1273 { |
|
1274 for (int i = 0; i < nr; i++) |
|
1275 work[i] = 0.; |
|
1276 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1277 work[b.ridx(i)] = b.data(i); |
|
1278 |
|
1279 for (int k = nr-1; k >= 0; k--) |
|
1280 { |
|
1281 int iidx = q_perm[k]; |
|
1282 if (work[iidx] != 0.) |
|
1283 { |
|
1284 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1285 { |
|
1286 err = -2; |
|
1287 goto triangular_error; |
|
1288 } |
|
1289 |
|
1290 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1291 work[iidx] = tmp; |
|
1292 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1293 { |
|
1294 int idx2 = q_perm[ridx(i)]; |
|
1295 work[idx2] = |
|
1296 work[idx2] - tmp * data(i); |
|
1297 } |
|
1298 } |
|
1299 } |
|
1300 |
|
1301 // Count non-zeros in work vector and adjust space in |
|
1302 // retval if needed |
|
1303 int new_nnz = 0; |
|
1304 for (int i = 0; i < nr; i++) |
|
1305 if (work[i] != 0.) |
|
1306 new_nnz++; |
|
1307 |
|
1308 if (ii + new_nnz > x_nz) |
|
1309 { |
|
1310 // Resize the sparse matrix |
|
1311 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1312 retval.change_capacity (sz); |
|
1313 x_nz = sz; |
|
1314 } |
|
1315 |
|
1316 for (int i = 0; i < nr; i++) |
|
1317 if (work[p_perm[i]] != 0.) |
|
1318 { |
|
1319 retval.xridx(ii) = i; |
|
1320 retval.xdata(ii++) = work[p_perm[i]]; |
|
1321 } |
|
1322 retval.xcidx(j+1) = ii; |
|
1323 } |
|
1324 |
|
1325 retval.maybe_compress (); |
|
1326 |
|
1327 // Calculation of 1-norm of inv(*this) |
|
1328 for (int i = 0; i < nr; i++) |
|
1329 work[i] = 0.; |
|
1330 |
|
1331 for (int j = 0; j < nr; j++) |
|
1332 { |
|
1333 work[q_perm[j]] = 1.; |
|
1334 |
|
1335 for (int k = j; k >= 0; k--) |
|
1336 { |
|
1337 int iidx = q_perm[k]; |
|
1338 |
|
1339 if (work[iidx] != 0.) |
|
1340 { |
|
1341 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1342 work[iidx] = tmp; |
|
1343 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1344 { |
|
1345 int idx2 = q_perm[ridx(i)]; |
|
1346 work[idx2] = work[idx2] - tmp * data(i); |
|
1347 } |
|
1348 } |
|
1349 } |
|
1350 double atmp = 0; |
|
1351 for (int i = 0; i < j+1; i++) |
|
1352 { |
|
1353 atmp += ::abs(work[i]); |
|
1354 work[i] = 0.; |
|
1355 } |
|
1356 if (atmp > ainvnorm) |
|
1357 ainvnorm = atmp; |
|
1358 } |
|
1359 } |
|
1360 else |
|
1361 { |
|
1362 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1363 |
|
1364 for (int j = 0; j < b_nc; j++) |
|
1365 { |
|
1366 for (int i = 0; i < nr; i++) |
|
1367 work[i] = 0.; |
|
1368 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1369 work[b.ridx(i)] = b.data(i); |
|
1370 |
|
1371 for (int k = nr-1; k >= 0; k--) |
|
1372 { |
|
1373 if (work[k] != 0.) |
|
1374 { |
|
1375 if (ridx(cidx(k+1)-1) != k) |
|
1376 { |
|
1377 err = -2; |
|
1378 goto triangular_error; |
|
1379 } |
|
1380 |
|
1381 Complex tmp = work[k] / data(cidx(k+1)-1); |
|
1382 work[k] = tmp; |
|
1383 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1384 { |
|
1385 int iidx = ridx(i); |
|
1386 work[iidx] = work[iidx] - tmp * data(i); |
|
1387 } |
|
1388 } |
|
1389 } |
|
1390 |
|
1391 // Count non-zeros in work vector and adjust space in |
|
1392 // retval if needed |
|
1393 int new_nnz = 0; |
|
1394 for (int i = 0; i < nr; i++) |
|
1395 if (work[i] != 0.) |
|
1396 new_nnz++; |
|
1397 |
|
1398 if (ii + new_nnz > x_nz) |
|
1399 { |
|
1400 // Resize the sparse matrix |
|
1401 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1402 retval.change_capacity (sz); |
|
1403 x_nz = sz; |
|
1404 } |
|
1405 |
|
1406 for (int i = 0; i < nr; i++) |
|
1407 if (work[i] != 0.) |
|
1408 { |
|
1409 retval.xridx(ii) = i; |
|
1410 retval.xdata(ii++) = work[i]; |
|
1411 } |
|
1412 retval.xcidx(j+1) = ii; |
|
1413 } |
|
1414 |
|
1415 retval.maybe_compress (); |
|
1416 |
|
1417 // Calculation of 1-norm of inv(*this) |
|
1418 for (int i = 0; i < nr; i++) |
|
1419 work[i] = 0.; |
|
1420 |
|
1421 for (int j = 0; j < nr; j++) |
|
1422 { |
|
1423 work[j] = 1.; |
|
1424 |
|
1425 for (int k = j; k >= 0; k--) |
|
1426 { |
|
1427 if (work[k] != 0.) |
|
1428 { |
|
1429 Complex tmp = work[k] / data(cidx(k+1)-1); |
|
1430 work[k] = tmp; |
|
1431 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1432 { |
|
1433 int iidx = ridx(i); |
|
1434 work[iidx] = work[iidx] - tmp * data(i); |
|
1435 } |
|
1436 } |
|
1437 } |
|
1438 double atmp = 0; |
|
1439 for (int i = 0; i < j+1; i++) |
|
1440 { |
|
1441 atmp += ::abs(work[i]); |
|
1442 work[i] = 0.; |
|
1443 } |
|
1444 if (atmp > ainvnorm) |
|
1445 ainvnorm = atmp; |
|
1446 } |
|
1447 } |
|
1448 |
|
1449 rcond = 1. / ainvnorm / anorm; |
|
1450 |
|
1451 triangular_error: |
|
1452 if (err != 0) |
|
1453 { |
|
1454 if (sing_handler) |
|
1455 sing_handler (rcond); |
|
1456 else |
|
1457 (*current_liboctave_error_handler) |
|
1458 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
1459 rcond); |
|
1460 } |
|
1461 |
|
1462 volatile double rcond_plus_one = rcond + 1.0; |
|
1463 |
|
1464 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
1465 { |
|
1466 err = -2; |
|
1467 |
|
1468 if (sing_handler) |
|
1469 sing_handler (rcond); |
|
1470 else |
|
1471 (*current_liboctave_error_handler) |
|
1472 ("matrix singular to machine precision, rcond = %g", |
|
1473 rcond); |
|
1474 } |
|
1475 } |
|
1476 else |
|
1477 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1478 } |
|
1479 return retval; |
|
1480 } |
|
1481 |
|
1482 ComplexMatrix |
|
1483 SparseComplexMatrix::utsolve (SparseType &mattype, const ComplexMatrix& b, |
|
1484 int& err, double& rcond, |
|
1485 solve_singularity_handler sing_handler) const |
|
1486 { |
|
1487 ComplexMatrix retval; |
|
1488 |
|
1489 int nr = rows (); |
|
1490 int nc = cols (); |
|
1491 err = 0; |
|
1492 |
|
1493 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1494 (*current_liboctave_error_handler) |
|
1495 ("matrix dimension mismatch solution of linear equations"); |
|
1496 else |
|
1497 { |
|
1498 // Print spparms("spumoni") info if requested |
|
1499 int typ = mattype.type (); |
|
1500 mattype.info (); |
|
1501 |
|
1502 if (typ == SparseType::Permuted_Upper || |
|
1503 typ == SparseType::Upper) |
|
1504 { |
|
1505 double anorm = 0.; |
|
1506 double ainvnorm = 0.; |
|
1507 int b_nc = b.cols (); |
|
1508 rcond = 0.; |
|
1509 |
|
1510 // Calculate the 1-norm of matrix for rcond calculation |
|
1511 for (int j = 0; j < nr; j++) |
|
1512 { |
|
1513 double atmp = 0.; |
|
1514 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1515 atmp += ::abs(data(i)); |
|
1516 if (atmp > anorm) |
|
1517 anorm = atmp; |
|
1518 } |
|
1519 |
|
1520 if (typ == SparseType::Permuted_Upper) |
|
1521 { |
|
1522 retval.resize (b.rows (), b.cols ()); |
|
1523 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1524 int *p_perm = mattype.triangular_row_perm (); |
|
1525 int *q_perm = mattype.triangular_col_perm (); |
|
1526 |
|
1527 (*current_liboctave_warning_handler) |
|
1528 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1529 |
|
1530 for (int j = 0; j < b_nc; j++) |
|
1531 { |
|
1532 for (int i = 0; i < nr; i++) |
|
1533 work[i] = b(i,j); |
|
1534 |
|
1535 for (int k = nr-1; k >= 0; k--) |
|
1536 { |
|
1537 int iidx = q_perm[k]; |
|
1538 if (work[iidx] != 0.) |
|
1539 { |
|
1540 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1541 { |
|
1542 err = -2; |
|
1543 goto triangular_error; |
|
1544 } |
|
1545 |
|
1546 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1547 work[iidx] = tmp; |
|
1548 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1549 { |
|
1550 int idx2 = q_perm[ridx(i)]; |
|
1551 work[idx2] = |
|
1552 work[idx2] - tmp * data(i); |
|
1553 } |
|
1554 } |
|
1555 } |
|
1556 |
|
1557 for (int i = 0; i < nr; i++) |
|
1558 retval (i, j) = work[p_perm[i]]; |
|
1559 |
|
1560 } |
|
1561 |
|
1562 // Calculation of 1-norm of inv(*this) |
|
1563 for (int i = 0; i < nr; i++) |
|
1564 work[i] = 0.; |
|
1565 |
|
1566 for (int j = 0; j < nr; j++) |
|
1567 { |
|
1568 work[q_perm[j]] = 1.; |
|
1569 |
|
1570 for (int k = j; k >= 0; k--) |
|
1571 { |
|
1572 int iidx = q_perm[k]; |
|
1573 |
|
1574 if (work[iidx] != 0.) |
|
1575 { |
|
1576 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1577 work[iidx] = tmp; |
|
1578 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1579 { |
|
1580 int idx2 = q_perm[ridx(i)]; |
|
1581 work[idx2] = work[idx2] - tmp * data(i); |
|
1582 } |
|
1583 } |
|
1584 } |
|
1585 double atmp = 0; |
|
1586 for (int i = 0; i < j+1; i++) |
|
1587 { |
|
1588 atmp += ::abs(work[i]); |
|
1589 work[i] = 0.; |
|
1590 } |
|
1591 if (atmp > ainvnorm) |
|
1592 ainvnorm = atmp; |
|
1593 } |
|
1594 } |
|
1595 else |
|
1596 { |
|
1597 retval = b; |
|
1598 Complex *x_vec = retval.fortran_vec (); |
|
1599 |
|
1600 for (int j = 0; j < b_nc; j++) |
|
1601 { |
|
1602 int offset = j * nr; |
|
1603 for (int k = nr-1; k >= 0; k--) |
|
1604 { |
|
1605 if (x_vec[k+offset] != 0.) |
|
1606 { |
|
1607 if (ridx(cidx(k+1)-1) != k) |
|
1608 { |
|
1609 err = -2; |
|
1610 goto triangular_error; |
|
1611 } |
|
1612 |
|
1613 Complex tmp = x_vec[k+offset] / |
|
1614 data(cidx(k+1)-1); |
|
1615 x_vec[k+offset] = tmp; |
|
1616 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1617 { |
|
1618 int iidx = ridx(i); |
|
1619 x_vec[iidx+offset] = |
|
1620 x_vec[iidx+offset] - tmp * data(i); |
|
1621 } |
|
1622 } |
|
1623 } |
|
1624 } |
|
1625 |
|
1626 // Calculation of 1-norm of inv(*this) |
|
1627 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1628 for (int i = 0; i < nr; i++) |
|
1629 work[i] = 0.; |
|
1630 |
|
1631 for (int j = 0; j < nr; j++) |
|
1632 { |
|
1633 work[j] = 1.; |
|
1634 |
|
1635 for (int k = j; k >= 0; k--) |
|
1636 { |
|
1637 if (work[k] != 0.) |
|
1638 { |
|
1639 Complex tmp = work[k] / data(cidx(k+1)-1); |
|
1640 work[k] = tmp; |
|
1641 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1642 { |
|
1643 int iidx = ridx(i); |
|
1644 work[iidx] = work[iidx] - tmp * data(i); |
|
1645 } |
|
1646 } |
|
1647 } |
|
1648 double atmp = 0; |
|
1649 for (int i = 0; i < j+1; i++) |
|
1650 { |
|
1651 atmp += ::abs(work[i]); |
|
1652 work[i] = 0.; |
|
1653 } |
|
1654 if (atmp > ainvnorm) |
|
1655 ainvnorm = atmp; |
|
1656 } |
|
1657 } |
|
1658 |
|
1659 rcond = 1. / ainvnorm / anorm; |
|
1660 |
|
1661 triangular_error: |
|
1662 if (err != 0) |
|
1663 { |
|
1664 if (sing_handler) |
|
1665 sing_handler (rcond); |
|
1666 else |
|
1667 (*current_liboctave_error_handler) |
|
1668 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
1669 rcond); |
|
1670 } |
|
1671 |
|
1672 volatile double rcond_plus_one = rcond + 1.0; |
|
1673 |
|
1674 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
1675 { |
|
1676 err = -2; |
|
1677 |
|
1678 if (sing_handler) |
|
1679 sing_handler (rcond); |
|
1680 else |
|
1681 (*current_liboctave_error_handler) |
|
1682 ("matrix singular to machine precision, rcond = %g", |
|
1683 rcond); |
|
1684 } |
|
1685 } |
|
1686 else |
|
1687 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1688 } |
|
1689 |
|
1690 return retval; |
|
1691 } |
|
1692 |
|
1693 SparseComplexMatrix |
|
1694 SparseComplexMatrix::utsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
1695 int& err, double& rcond, |
|
1696 solve_singularity_handler sing_handler) const |
|
1697 { |
|
1698 SparseComplexMatrix retval; |
|
1699 |
|
1700 int nr = rows (); |
|
1701 int nc = cols (); |
|
1702 err = 0; |
|
1703 |
|
1704 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1705 (*current_liboctave_error_handler) |
|
1706 ("matrix dimension mismatch solution of linear equations"); |
|
1707 else |
|
1708 { |
|
1709 // Print spparms("spumoni") info if requested |
|
1710 int typ = mattype.type (); |
|
1711 mattype.info (); |
|
1712 |
|
1713 if (typ == SparseType::Permuted_Upper || |
|
1714 typ == SparseType::Upper) |
|
1715 { |
|
1716 double anorm = 0.; |
|
1717 double ainvnorm = 0.; |
|
1718 rcond = 0.; |
|
1719 |
|
1720 // Calculate the 1-norm of matrix for rcond calculation |
|
1721 for (int j = 0; j < nr; j++) |
|
1722 { |
|
1723 double atmp = 0.; |
|
1724 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1725 atmp += ::abs(data(i)); |
|
1726 if (atmp > anorm) |
|
1727 anorm = atmp; |
|
1728 } |
|
1729 |
|
1730 int b_nr = b.rows (); |
|
1731 int b_nc = b.cols (); |
|
1732 int b_nz = b.nnz (); |
|
1733 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
1734 retval.xcidx(0) = 0; |
|
1735 int ii = 0; |
|
1736 int x_nz = b_nz; |
|
1737 |
|
1738 if (typ == SparseType::Permuted_Upper) |
|
1739 { |
|
1740 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1741 int *p_perm = mattype.triangular_row_perm (); |
|
1742 int *q_perm = mattype.triangular_col_perm (); |
|
1743 |
|
1744 (*current_liboctave_warning_handler) |
|
1745 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1746 |
|
1747 for (int j = 0; j < b_nc; j++) |
|
1748 { |
|
1749 for (int i = 0; i < nr; i++) |
|
1750 work[i] = 0.; |
|
1751 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1752 work[b.ridx(i)] = b.data(i); |
|
1753 |
|
1754 for (int k = nr-1; k >= 0; k--) |
|
1755 { |
|
1756 int iidx = q_perm[k]; |
|
1757 if (work[iidx] != 0.) |
|
1758 { |
|
1759 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1760 { |
|
1761 err = -2; |
|
1762 goto triangular_error; |
|
1763 } |
|
1764 |
|
1765 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1766 work[iidx] = tmp; |
|
1767 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1768 { |
|
1769 int idx2 = q_perm[ridx(i)]; |
|
1770 work[idx2] = |
|
1771 work[idx2] - tmp * data(i); |
|
1772 } |
|
1773 } |
|
1774 } |
|
1775 |
|
1776 // Count non-zeros in work vector and adjust space in |
|
1777 // retval if needed |
|
1778 int new_nnz = 0; |
|
1779 for (int i = 0; i < nr; i++) |
|
1780 if (work[i] != 0.) |
|
1781 new_nnz++; |
|
1782 |
|
1783 if (ii + new_nnz > x_nz) |
|
1784 { |
|
1785 // Resize the sparse matrix |
|
1786 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1787 retval.change_capacity (sz); |
|
1788 x_nz = sz; |
|
1789 } |
|
1790 |
|
1791 for (int i = 0; i < nr; i++) |
|
1792 if (work[p_perm[i]] != 0.) |
|
1793 { |
|
1794 retval.xridx(ii) = i; |
|
1795 retval.xdata(ii++) = work[p_perm[i]]; |
|
1796 } |
|
1797 retval.xcidx(j+1) = ii; |
|
1798 } |
|
1799 |
|
1800 retval.maybe_compress (); |
|
1801 |
|
1802 // Calculation of 1-norm of inv(*this) |
|
1803 for (int i = 0; i < nr; i++) |
|
1804 work[i] = 0.; |
|
1805 |
|
1806 for (int j = 0; j < nr; j++) |
|
1807 { |
|
1808 work[q_perm[j]] = 1.; |
|
1809 |
|
1810 for (int k = j; k >= 0; k--) |
|
1811 { |
|
1812 int iidx = q_perm[k]; |
|
1813 |
|
1814 if (work[iidx] != 0.) |
|
1815 { |
|
1816 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1817 work[iidx] = tmp; |
|
1818 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1819 { |
|
1820 int idx2 = q_perm[ridx(i)]; |
|
1821 work[idx2] = work[idx2] - tmp * data(i); |
|
1822 } |
|
1823 } |
|
1824 } |
|
1825 double atmp = 0; |
|
1826 for (int i = 0; i < j+1; i++) |
|
1827 { |
|
1828 atmp += ::abs(work[i]); |
|
1829 work[i] = 0.; |
|
1830 } |
|
1831 if (atmp > ainvnorm) |
|
1832 ainvnorm = atmp; |
|
1833 } |
|
1834 } |
|
1835 else |
|
1836 { |
|
1837 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1838 |
|
1839 for (int j = 0; j < b_nc; j++) |
|
1840 { |
|
1841 for (int i = 0; i < nr; i++) |
|
1842 work[i] = 0.; |
|
1843 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1844 work[b.ridx(i)] = b.data(i); |
|
1845 |
|
1846 for (int k = nr-1; k >= 0; k--) |
|
1847 { |
|
1848 if (work[k] != 0.) |
|
1849 { |
|
1850 if (ridx(cidx(k+1)-1) != k) |
|
1851 { |
|
1852 err = -2; |
|
1853 goto triangular_error; |
|
1854 } |
|
1855 |
|
1856 Complex tmp = work[k] / data(cidx(k+1)-1); |
|
1857 work[k] = tmp; |
|
1858 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1859 { |
|
1860 int iidx = ridx(i); |
|
1861 work[iidx] = work[iidx] - tmp * data(i); |
|
1862 } |
|
1863 } |
|
1864 } |
|
1865 |
|
1866 // Count non-zeros in work vector and adjust space in |
|
1867 // retval if needed |
|
1868 int new_nnz = 0; |
|
1869 for (int i = 0; i < nr; i++) |
|
1870 if (work[i] != 0.) |
|
1871 new_nnz++; |
|
1872 |
|
1873 if (ii + new_nnz > x_nz) |
|
1874 { |
|
1875 // Resize the sparse matrix |
|
1876 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1877 retval.change_capacity (sz); |
|
1878 x_nz = sz; |
|
1879 } |
|
1880 |
|
1881 for (int i = 0; i < nr; i++) |
|
1882 if (work[i] != 0.) |
|
1883 { |
|
1884 retval.xridx(ii) = i; |
|
1885 retval.xdata(ii++) = work[i]; |
|
1886 } |
|
1887 retval.xcidx(j+1) = ii; |
|
1888 } |
|
1889 |
|
1890 retval.maybe_compress (); |
|
1891 |
|
1892 // Calculation of 1-norm of inv(*this) |
|
1893 for (int i = 0; i < nr; i++) |
|
1894 work[i] = 0.; |
|
1895 |
|
1896 for (int j = 0; j < nr; j++) |
|
1897 { |
|
1898 work[j] = 1.; |
|
1899 |
|
1900 for (int k = j; k >= 0; k--) |
|
1901 { |
|
1902 if (work[k] != 0.) |
|
1903 { |
|
1904 Complex tmp = work[k] / data(cidx(k+1)-1); |
|
1905 work[k] = tmp; |
|
1906 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1907 { |
|
1908 int iidx = ridx(i); |
|
1909 work[iidx] = work[iidx] - tmp * data(i); |
|
1910 } |
|
1911 } |
|
1912 } |
|
1913 double atmp = 0; |
|
1914 for (int i = 0; i < j+1; i++) |
|
1915 { |
|
1916 atmp += ::abs(work[i]); |
|
1917 work[i] = 0.; |
|
1918 } |
|
1919 if (atmp > ainvnorm) |
|
1920 ainvnorm = atmp; |
|
1921 } |
|
1922 } |
|
1923 |
|
1924 rcond = 1. / ainvnorm / anorm; |
|
1925 |
|
1926 triangular_error: |
|
1927 if (err != 0) |
|
1928 { |
|
1929 if (sing_handler) |
|
1930 sing_handler (rcond); |
|
1931 else |
|
1932 (*current_liboctave_error_handler) |
|
1933 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
1934 rcond); |
|
1935 } |
|
1936 |
|
1937 volatile double rcond_plus_one = rcond + 1.0; |
|
1938 |
|
1939 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
1940 { |
|
1941 err = -2; |
|
1942 |
|
1943 if (sing_handler) |
|
1944 sing_handler (rcond); |
|
1945 else |
|
1946 (*current_liboctave_error_handler) |
|
1947 ("matrix singular to machine precision, rcond = %g", |
|
1948 rcond); |
|
1949 } |
|
1950 } |
|
1951 else |
|
1952 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1953 } |
|
1954 |
|
1955 return retval; |
|
1956 } |
|
1957 |
|
1958 ComplexMatrix |
|
1959 SparseComplexMatrix::ltsolve (SparseType &mattype, const Matrix& b, int& err, |
|
1960 double& rcond, solve_singularity_handler sing_handler) const |
|
1961 { |
|
1962 ComplexMatrix retval; |
|
1963 |
|
1964 int nr = rows (); |
|
1965 int nc = cols (); |
|
1966 err = 0; |
|
1967 |
|
1968 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1969 (*current_liboctave_error_handler) |
|
1970 ("matrix dimension mismatch solution of linear equations"); |
|
1971 else |
|
1972 { |
|
1973 // Print spparms("spumoni") info if requested |
|
1974 int typ = mattype.type (); |
|
1975 mattype.info (); |
|
1976 |
|
1977 if (typ == SparseType::Permuted_Lower || |
|
1978 typ == SparseType::Lower) |
|
1979 { |
|
1980 double anorm = 0.; |
|
1981 double ainvnorm = 0.; |
|
1982 int b_cols = b.cols (); |
|
1983 rcond = 0.; |
|
1984 |
|
1985 // Calculate the 1-norm of matrix for rcond calculation |
|
1986 for (int j = 0; j < nr; j++) |
|
1987 { |
|
1988 double atmp = 0.; |
|
1989 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1990 atmp += ::abs(data(i)); |
|
1991 if (atmp > anorm) |
|
1992 anorm = atmp; |
|
1993 } |
|
1994 |
|
1995 if (typ == SparseType::Permuted_Lower) |
|
1996 { |
|
1997 retval.resize (b.rows (), b.cols ()); |
|
1998 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1999 int *p_perm = mattype.triangular_row_perm (); |
|
2000 int *q_perm = mattype.triangular_col_perm (); |
|
2001 |
|
2002 (*current_liboctave_warning_handler) |
|
2003 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2004 |
|
2005 for (int j = 0; j < b_cols; j++) |
|
2006 { |
|
2007 for (int i = 0; i < nr; i++) |
|
2008 work[i] = b(i,j); |
|
2009 |
|
2010 for (int k = 0; k < nr; k++) |
|
2011 { |
|
2012 int iidx = q_perm[k]; |
|
2013 if (work[iidx] != 0.) |
|
2014 { |
|
2015 if (ridx(cidx(iidx)) != iidx) |
|
2016 { |
|
2017 err = -2; |
|
2018 goto triangular_error; |
|
2019 } |
|
2020 |
|
2021 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2022 work[iidx] = tmp; |
|
2023 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2024 { |
|
2025 int idx2 = q_perm[ridx(i)]; |
|
2026 work[idx2] = |
|
2027 work[idx2] - tmp * data(i); |
|
2028 } |
|
2029 } |
|
2030 } |
|
2031 |
|
2032 for (int i = 0; i < nr; i++) |
|
2033 retval (i, j) = work[p_perm[i]]; |
|
2034 |
|
2035 } |
|
2036 |
|
2037 // Calculation of 1-norm of inv(*this) |
|
2038 for (int i = 0; i < nr; i++) |
|
2039 work[i] = 0.; |
|
2040 |
|
2041 for (int j = 0; j < nr; j++) |
|
2042 { |
|
2043 work[q_perm[j]] = 1.; |
|
2044 |
|
2045 for (int k = 0; k < nr; k++) |
|
2046 { |
|
2047 int iidx = q_perm[k]; |
|
2048 |
|
2049 if (work[iidx] != 0.) |
|
2050 { |
|
2051 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2052 work[iidx] = tmp; |
|
2053 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2054 { |
|
2055 int idx2 = q_perm[ridx(i)]; |
|
2056 work[idx2] = work[idx2] - tmp * data(i); |
|
2057 } |
|
2058 } |
|
2059 } |
|
2060 double atmp = 0; |
|
2061 for (int i = 0; i < j+1; i++) |
|
2062 { |
|
2063 atmp += ::abs(work[i]); |
|
2064 work[i] = 0.; |
|
2065 } |
|
2066 if (atmp > ainvnorm) |
|
2067 ainvnorm = atmp; |
|
2068 } |
|
2069 } |
|
2070 else |
|
2071 { |
|
2072 retval = ComplexMatrix (b); |
|
2073 Complex *x_vec = retval.fortran_vec (); |
|
2074 |
|
2075 for (int j = 0; j < b_cols; j++) |
|
2076 { |
|
2077 int offset = j * nr; |
|
2078 for (int k = 0; k < nr; k++) |
|
2079 { |
|
2080 if (x_vec[k+offset] != 0.) |
|
2081 { |
|
2082 if (ridx(cidx(k)) != k) |
|
2083 { |
|
2084 err = -2; |
|
2085 goto triangular_error; |
|
2086 } |
|
2087 |
|
2088 Complex tmp = x_vec[k+offset] / |
|
2089 data(cidx(k)); |
|
2090 x_vec[k+offset] = tmp; |
|
2091 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2092 { |
|
2093 int iidx = ridx(i); |
|
2094 x_vec[iidx+offset] = |
|
2095 x_vec[iidx+offset] - tmp * data(i); |
|
2096 } |
|
2097 } |
|
2098 } |
|
2099 } |
|
2100 |
|
2101 // Calculation of 1-norm of inv(*this) |
|
2102 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2103 for (int i = 0; i < nr; i++) |
|
2104 work[i] = 0.; |
|
2105 |
|
2106 for (int j = 0; j < nr; j++) |
|
2107 { |
|
2108 work[j] = 1.; |
|
2109 |
|
2110 for (int k = j; k < nr; k++) |
|
2111 { |
|
2112 |
|
2113 if (work[k] != 0.) |
|
2114 { |
|
2115 Complex tmp = work[k] / data(cidx(k)); |
|
2116 work[k] = tmp; |
|
2117 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2118 { |
|
2119 int iidx = ridx(i); |
|
2120 work[iidx] = work[iidx] - tmp * data(i); |
|
2121 } |
|
2122 } |
|
2123 } |
|
2124 double atmp = 0; |
|
2125 for (int i = j; i < nr; i++) |
|
2126 { |
|
2127 atmp += ::abs(work[i]); |
|
2128 work[i] = 0.; |
|
2129 } |
|
2130 if (atmp > ainvnorm) |
|
2131 ainvnorm = atmp; |
|
2132 } |
|
2133 } |
|
2134 |
|
2135 rcond = 1. / ainvnorm / anorm; |
|
2136 |
|
2137 triangular_error: |
|
2138 if (err != 0) |
|
2139 { |
|
2140 if (sing_handler) |
|
2141 sing_handler (rcond); |
|
2142 else |
|
2143 (*current_liboctave_error_handler) |
|
2144 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2145 rcond); |
|
2146 } |
|
2147 |
|
2148 volatile double rcond_plus_one = rcond + 1.0; |
|
2149 |
|
2150 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2151 { |
|
2152 err = -2; |
|
2153 |
|
2154 if (sing_handler) |
|
2155 sing_handler (rcond); |
|
2156 else |
|
2157 (*current_liboctave_error_handler) |
|
2158 ("matrix singular to machine precision, rcond = %g", |
|
2159 rcond); |
|
2160 } |
|
2161 } |
|
2162 else |
|
2163 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2164 } |
|
2165 |
|
2166 return retval; |
|
2167 } |
|
2168 |
|
2169 SparseComplexMatrix |
|
2170 SparseComplexMatrix::ltsolve (SparseType &mattype, const SparseMatrix& b, |
|
2171 int& err, double& rcond, |
|
2172 solve_singularity_handler sing_handler) const |
|
2173 { |
|
2174 SparseComplexMatrix retval; |
|
2175 |
|
2176 int nr = rows (); |
|
2177 int nc = cols (); |
|
2178 err = 0; |
|
2179 |
|
2180 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2181 (*current_liboctave_error_handler) |
|
2182 ("matrix dimension mismatch solution of linear equations"); |
|
2183 else |
|
2184 { |
|
2185 // Print spparms("spumoni") info if requested |
|
2186 int typ = mattype.type (); |
|
2187 mattype.info (); |
|
2188 |
|
2189 if (typ == SparseType::Permuted_Lower || |
|
2190 typ == SparseType::Lower) |
|
2191 { |
|
2192 double anorm = 0.; |
|
2193 double ainvnorm = 0.; |
|
2194 rcond = 0.; |
|
2195 |
|
2196 // Calculate the 1-norm of matrix for rcond calculation |
|
2197 for (int j = 0; j < nr; j++) |
|
2198 { |
|
2199 double atmp = 0.; |
|
2200 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2201 atmp += ::abs(data(i)); |
|
2202 if (atmp > anorm) |
|
2203 anorm = atmp; |
|
2204 } |
|
2205 |
|
2206 int b_nr = b.rows (); |
|
2207 int b_nc = b.cols (); |
|
2208 int b_nz = b.nnz (); |
|
2209 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
2210 retval.xcidx(0) = 0; |
|
2211 int ii = 0; |
|
2212 int x_nz = b_nz; |
|
2213 |
|
2214 if (typ == SparseType::Permuted_Lower) |
|
2215 { |
|
2216 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2217 int *p_perm = mattype.triangular_row_perm (); |
|
2218 int *q_perm = mattype.triangular_col_perm (); |
|
2219 |
|
2220 (*current_liboctave_warning_handler) |
|
2221 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2222 |
|
2223 for (int j = 0; j < b_nc; j++) |
|
2224 { |
|
2225 for (int i = 0; i < nr; i++) |
|
2226 work[i] = 0.; |
|
2227 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2228 work[b.ridx(i)] = b.data(i); |
|
2229 |
|
2230 for (int k = 0; k < nr; k++) |
|
2231 { |
|
2232 int iidx = q_perm[k]; |
|
2233 if (work[iidx] != 0.) |
|
2234 { |
|
2235 if (ridx(cidx(iidx)) != iidx) |
|
2236 { |
|
2237 err = -2; |
|
2238 goto triangular_error; |
|
2239 } |
|
2240 |
|
2241 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2242 work[iidx] = tmp; |
|
2243 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2244 { |
|
2245 int idx2 = q_perm[ridx(i)]; |
|
2246 work[idx2] = |
|
2247 work[idx2] - tmp * data(i); |
|
2248 } |
|
2249 } |
|
2250 } |
|
2251 |
|
2252 // Count non-zeros in work vector and adjust space in |
|
2253 // retval if needed |
|
2254 int new_nnz = 0; |
|
2255 for (int i = 0; i < nr; i++) |
|
2256 if (work[i] != 0.) |
|
2257 new_nnz++; |
|
2258 |
|
2259 if (ii + new_nnz > x_nz) |
|
2260 { |
|
2261 // Resize the sparse matrix |
|
2262 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2263 retval.change_capacity (sz); |
|
2264 x_nz = sz; |
|
2265 } |
|
2266 |
|
2267 for (int i = 0; i < nr; i++) |
|
2268 if (work[p_perm[i]] != 0.) |
|
2269 { |
|
2270 retval.xridx(ii) = i; |
|
2271 retval.xdata(ii++) = work[p_perm[i]]; |
|
2272 } |
|
2273 retval.xcidx(j+1) = ii; |
|
2274 } |
|
2275 |
|
2276 retval.maybe_compress (); |
|
2277 |
|
2278 // Calculation of 1-norm of inv(*this) |
|
2279 for (int i = 0; i < nr; i++) |
|
2280 work[i] = 0.; |
|
2281 |
|
2282 for (int j = 0; j < nr; j++) |
|
2283 { |
|
2284 work[q_perm[j]] = 1.; |
|
2285 |
|
2286 for (int k = 0; k < nr; k++) |
|
2287 { |
|
2288 int iidx = q_perm[k]; |
|
2289 |
|
2290 if (work[iidx] != 0.) |
|
2291 { |
|
2292 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2293 work[iidx] = tmp; |
|
2294 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2295 { |
|
2296 int idx2 = q_perm[ridx(i)]; |
|
2297 work[idx2] = work[idx2] - tmp * data(i); |
|
2298 } |
|
2299 } |
|
2300 } |
|
2301 double atmp = 0; |
|
2302 for (int i = 0; i < j+1; i++) |
|
2303 { |
|
2304 atmp += ::abs(work[i]); |
|
2305 work[i] = 0.; |
|
2306 } |
|
2307 if (atmp > ainvnorm) |
|
2308 ainvnorm = atmp; |
|
2309 } |
|
2310 } |
|
2311 else |
|
2312 { |
|
2313 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2314 |
|
2315 for (int j = 0; j < b_nc; j++) |
|
2316 { |
|
2317 for (int i = 0; i < nr; i++) |
|
2318 work[i] = 0.; |
|
2319 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2320 work[b.ridx(i)] = b.data(i); |
|
2321 |
|
2322 for (int k = 0; k < nr; k++) |
|
2323 { |
|
2324 if (work[k] != 0.) |
|
2325 { |
|
2326 if (ridx(cidx(k)) != k) |
|
2327 { |
|
2328 err = -2; |
|
2329 goto triangular_error; |
|
2330 } |
|
2331 |
|
2332 Complex tmp = work[k] / data(cidx(k)); |
|
2333 work[k] = tmp; |
|
2334 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2335 { |
|
2336 int iidx = ridx(i); |
|
2337 work[iidx] = work[iidx] - tmp * data(i); |
|
2338 } |
|
2339 } |
|
2340 } |
|
2341 |
|
2342 // Count non-zeros in work vector and adjust space in |
|
2343 // retval if needed |
|
2344 int new_nnz = 0; |
|
2345 for (int i = 0; i < nr; i++) |
|
2346 if (work[i] != 0.) |
|
2347 new_nnz++; |
|
2348 |
|
2349 if (ii + new_nnz > x_nz) |
|
2350 { |
|
2351 // Resize the sparse matrix |
|
2352 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2353 retval.change_capacity (sz); |
|
2354 x_nz = sz; |
|
2355 } |
|
2356 |
|
2357 for (int i = 0; i < nr; i++) |
|
2358 if (work[i] != 0.) |
|
2359 { |
|
2360 retval.xridx(ii) = i; |
|
2361 retval.xdata(ii++) = work[i]; |
|
2362 } |
|
2363 retval.xcidx(j+1) = ii; |
|
2364 } |
|
2365 |
|
2366 retval.maybe_compress (); |
|
2367 |
|
2368 // Calculation of 1-norm of inv(*this) |
|
2369 for (int i = 0; i < nr; i++) |
|
2370 work[i] = 0.; |
|
2371 |
|
2372 for (int j = 0; j < nr; j++) |
|
2373 { |
|
2374 work[j] = 1.; |
|
2375 |
|
2376 for (int k = j; k < nr; k++) |
|
2377 { |
|
2378 |
|
2379 if (work[k] != 0.) |
|
2380 { |
|
2381 Complex tmp = work[k] / data(cidx(k)); |
|
2382 work[k] = tmp; |
|
2383 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2384 { |
|
2385 int iidx = ridx(i); |
|
2386 work[iidx] = work[iidx] - tmp * data(i); |
|
2387 } |
|
2388 } |
|
2389 } |
|
2390 double atmp = 0; |
|
2391 for (int i = j; i < nr; i++) |
|
2392 { |
|
2393 atmp += ::abs(work[i]); |
|
2394 work[i] = 0.; |
|
2395 } |
|
2396 if (atmp > ainvnorm) |
|
2397 ainvnorm = atmp; |
|
2398 } |
|
2399 |
|
2400 } |
|
2401 |
|
2402 rcond = 1. / ainvnorm / anorm; |
|
2403 |
|
2404 triangular_error: |
|
2405 if (err != 0) |
|
2406 { |
|
2407 if (sing_handler) |
|
2408 sing_handler (rcond); |
|
2409 else |
|
2410 (*current_liboctave_error_handler) |
|
2411 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2412 rcond); |
|
2413 } |
|
2414 |
|
2415 volatile double rcond_plus_one = rcond + 1.0; |
|
2416 |
|
2417 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2418 { |
|
2419 err = -2; |
|
2420 |
|
2421 if (sing_handler) |
|
2422 sing_handler (rcond); |
|
2423 else |
|
2424 (*current_liboctave_error_handler) |
|
2425 ("matrix singular to machine precision, rcond = %g", |
|
2426 rcond); |
|
2427 } |
|
2428 } |
|
2429 else |
|
2430 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2431 } |
|
2432 |
|
2433 return retval; |
|
2434 } |
|
2435 |
|
2436 ComplexMatrix |
|
2437 SparseComplexMatrix::ltsolve (SparseType &mattype, const ComplexMatrix& b, |
|
2438 int& err, double& rcond, |
|
2439 solve_singularity_handler sing_handler) const |
|
2440 { |
|
2441 ComplexMatrix retval; |
|
2442 |
|
2443 int nr = rows (); |
|
2444 int nc = cols (); |
|
2445 err = 0; |
|
2446 |
|
2447 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2448 (*current_liboctave_error_handler) |
|
2449 ("matrix dimension mismatch solution of linear equations"); |
|
2450 else |
|
2451 { |
|
2452 // Print spparms("spumoni") info if requested |
|
2453 int typ = mattype.type (); |
|
2454 mattype.info (); |
|
2455 |
|
2456 if (typ == SparseType::Permuted_Lower || |
|
2457 typ == SparseType::Lower) |
|
2458 { |
|
2459 double anorm = 0.; |
|
2460 double ainvnorm = 0.; |
|
2461 int b_nc = b.cols (); |
|
2462 rcond = 0.; |
|
2463 |
|
2464 // Calculate the 1-norm of matrix for rcond calculation |
|
2465 for (int j = 0; j < nr; j++) |
|
2466 { |
|
2467 double atmp = 0.; |
|
2468 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2469 atmp += ::abs(data(i)); |
|
2470 if (atmp > anorm) |
|
2471 anorm = atmp; |
|
2472 } |
|
2473 |
|
2474 if (typ == SparseType::Permuted_Lower) |
|
2475 { |
|
2476 retval.resize (b.rows (), b.cols ()); |
|
2477 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2478 int *p_perm = mattype.triangular_row_perm (); |
|
2479 int *q_perm = mattype.triangular_col_perm (); |
|
2480 |
|
2481 (*current_liboctave_warning_handler) |
|
2482 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2483 |
|
2484 for (int j = 0; j < b_nc; j++) |
|
2485 { |
|
2486 for (int i = 0; i < nr; i++) |
|
2487 work[i] = b(i,j); |
|
2488 |
|
2489 for (int k = 0; k < nr; k++) |
|
2490 { |
|
2491 int iidx = q_perm[k]; |
|
2492 if (work[iidx] != 0.) |
|
2493 { |
|
2494 if (ridx(cidx(iidx)) != iidx) |
|
2495 { |
|
2496 err = -2; |
|
2497 goto triangular_error; |
|
2498 } |
|
2499 |
|
2500 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2501 work[iidx] = tmp; |
|
2502 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2503 { |
|
2504 int idx2 = q_perm[ridx(i)]; |
|
2505 work[idx2] = |
|
2506 work[idx2] - tmp * data(i); |
|
2507 } |
|
2508 } |
|
2509 } |
|
2510 |
|
2511 for (int i = 0; i < nr; i++) |
|
2512 retval (i, j) = work[p_perm[i]]; |
|
2513 |
|
2514 } |
|
2515 |
|
2516 // Calculation of 1-norm of inv(*this) |
|
2517 for (int i = 0; i < nr; i++) |
|
2518 work[i] = 0.; |
|
2519 |
|
2520 for (int j = 0; j < nr; j++) |
|
2521 { |
|
2522 work[q_perm[j]] = 1.; |
|
2523 |
|
2524 for (int k = 0; k < nr; k++) |
|
2525 { |
|
2526 int iidx = q_perm[k]; |
|
2527 |
|
2528 if (work[iidx] != 0.) |
|
2529 { |
|
2530 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2531 work[iidx] = tmp; |
|
2532 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2533 { |
|
2534 int idx2 = q_perm[ridx(i)]; |
|
2535 work[idx2] = work[idx2] - tmp * data(i); |
|
2536 } |
|
2537 } |
|
2538 } |
|
2539 double atmp = 0; |
|
2540 for (int i = 0; i < j+1; i++) |
|
2541 { |
|
2542 atmp += ::abs(work[i]); |
|
2543 work[i] = 0.; |
|
2544 } |
|
2545 if (atmp > ainvnorm) |
|
2546 ainvnorm = atmp; |
|
2547 } |
|
2548 } |
|
2549 else |
|
2550 { |
|
2551 retval = b; |
|
2552 Complex *x_vec = retval.fortran_vec (); |
|
2553 |
|
2554 for (int j = 0; j < b_nc; j++) |
|
2555 { |
|
2556 int offset = j * nr; |
|
2557 for (int k = 0; k < nr; k++) |
|
2558 { |
|
2559 if (x_vec[k+offset] != 0.) |
|
2560 { |
|
2561 if (ridx(cidx(k)) != k) |
|
2562 { |
|
2563 err = -2; |
|
2564 goto triangular_error; |
|
2565 } |
|
2566 |
|
2567 Complex tmp = x_vec[k+offset] / |
|
2568 data(cidx(k)); |
|
2569 x_vec[k+offset] = tmp; |
|
2570 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2571 { |
|
2572 int iidx = ridx(i); |
|
2573 x_vec[iidx+offset] = |
|
2574 x_vec[iidx+offset] - tmp * data(i); |
|
2575 } |
|
2576 } |
|
2577 } |
|
2578 } |
|
2579 |
|
2580 // Calculation of 1-norm of inv(*this) |
|
2581 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2582 for (int i = 0; i < nr; i++) |
|
2583 work[i] = 0.; |
|
2584 |
|
2585 for (int j = 0; j < nr; j++) |
|
2586 { |
|
2587 work[j] = 1.; |
|
2588 |
|
2589 for (int k = j; k < nr; k++) |
|
2590 { |
|
2591 |
|
2592 if (work[k] != 0.) |
|
2593 { |
|
2594 Complex tmp = work[k] / data(cidx(k)); |
|
2595 work[k] = tmp; |
|
2596 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2597 { |
|
2598 int iidx = ridx(i); |
|
2599 work[iidx] = work[iidx] - tmp * data(i); |
|
2600 } |
|
2601 } |
|
2602 } |
|
2603 double atmp = 0; |
|
2604 for (int i = j; i < nr; i++) |
|
2605 { |
|
2606 atmp += ::abs(work[i]); |
|
2607 work[i] = 0.; |
|
2608 } |
|
2609 if (atmp > ainvnorm) |
|
2610 ainvnorm = atmp; |
|
2611 } |
|
2612 |
|
2613 } |
|
2614 |
|
2615 rcond = 1. / ainvnorm / anorm; |
|
2616 |
|
2617 triangular_error: |
|
2618 if (err != 0) |
|
2619 { |
|
2620 if (sing_handler) |
|
2621 sing_handler (rcond); |
|
2622 else |
|
2623 (*current_liboctave_error_handler) |
|
2624 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2625 rcond); |
|
2626 } |
|
2627 |
|
2628 volatile double rcond_plus_one = rcond + 1.0; |
|
2629 |
|
2630 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2631 { |
|
2632 err = -2; |
|
2633 |
|
2634 if (sing_handler) |
|
2635 sing_handler (rcond); |
|
2636 else |
|
2637 (*current_liboctave_error_handler) |
|
2638 ("matrix singular to machine precision, rcond = %g", |
|
2639 rcond); |
|
2640 } |
|
2641 } |
|
2642 else |
|
2643 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2644 } |
|
2645 |
|
2646 return retval; |
|
2647 } |
|
2648 |
|
2649 SparseComplexMatrix |
|
2650 SparseComplexMatrix::ltsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
2651 int& err, double& rcond, |
|
2652 solve_singularity_handler sing_handler) const |
|
2653 { |
|
2654 SparseComplexMatrix retval; |
|
2655 |
|
2656 int nr = rows (); |
|
2657 int nc = cols (); |
|
2658 err = 0; |
|
2659 |
|
2660 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2661 (*current_liboctave_error_handler) |
|
2662 ("matrix dimension mismatch solution of linear equations"); |
|
2663 else |
|
2664 { |
|
2665 // Print spparms("spumoni") info if requested |
|
2666 int typ = mattype.type (); |
|
2667 mattype.info (); |
|
2668 |
|
2669 if (typ == SparseType::Permuted_Lower || |
|
2670 typ == SparseType::Lower) |
|
2671 { |
|
2672 double anorm = 0.; |
|
2673 double ainvnorm = 0.; |
|
2674 rcond = 0.; |
|
2675 |
|
2676 // Calculate the 1-norm of matrix for rcond calculation |
|
2677 for (int j = 0; j < nr; j++) |
|
2678 { |
|
2679 double atmp = 0.; |
|
2680 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2681 atmp += ::abs(data(i)); |
|
2682 if (atmp > anorm) |
|
2683 anorm = atmp; |
|
2684 } |
|
2685 |
|
2686 int b_nr = b.rows (); |
|
2687 int b_nc = b.cols (); |
|
2688 int b_nz = b.nnz (); |
|
2689 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
2690 retval.xcidx(0) = 0; |
|
2691 int ii = 0; |
|
2692 int x_nz = b_nz; |
|
2693 |
|
2694 if (typ == SparseType::Permuted_Lower) |
|
2695 { |
|
2696 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2697 int *p_perm = mattype.triangular_row_perm (); |
|
2698 int *q_perm = mattype.triangular_col_perm (); |
|
2699 |
|
2700 (*current_liboctave_warning_handler) |
|
2701 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2702 |
|
2703 for (int j = 0; j < b_nc; j++) |
|
2704 { |
|
2705 for (int i = 0; i < nr; i++) |
|
2706 work[i] = 0.; |
|
2707 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2708 work[b.ridx(i)] = b.data(i); |
|
2709 |
|
2710 for (int k = 0; k < nr; k++) |
|
2711 { |
|
2712 int iidx = q_perm[k]; |
|
2713 if (work[iidx] != 0.) |
|
2714 { |
|
2715 if (ridx(cidx(iidx)) != iidx) |
|
2716 { |
|
2717 err = -2; |
|
2718 goto triangular_error; |
|
2719 } |
|
2720 |
|
2721 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2722 work[iidx] = tmp; |
|
2723 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2724 { |
|
2725 int idx2 = q_perm[ridx(i)]; |
|
2726 work[idx2] = |
|
2727 work[idx2] - tmp * data(i); |
|
2728 } |
|
2729 } |
|
2730 } |
|
2731 |
|
2732 // Count non-zeros in work vector and adjust space in |
|
2733 // retval if needed |
|
2734 int new_nnz = 0; |
|
2735 for (int i = 0; i < nr; i++) |
|
2736 if (work[i] != 0.) |
|
2737 new_nnz++; |
|
2738 |
|
2739 if (ii + new_nnz > x_nz) |
|
2740 { |
|
2741 // Resize the sparse matrix |
|
2742 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2743 retval.change_capacity (sz); |
|
2744 x_nz = sz; |
|
2745 } |
|
2746 |
|
2747 for (int i = 0; i < nr; i++) |
|
2748 if (work[p_perm[i]] != 0.) |
|
2749 { |
|
2750 retval.xridx(ii) = i; |
|
2751 retval.xdata(ii++) = work[p_perm[i]]; |
|
2752 } |
|
2753 retval.xcidx(j+1) = ii; |
|
2754 } |
|
2755 |
|
2756 retval.maybe_compress (); |
|
2757 |
|
2758 // Calculation of 1-norm of inv(*this) |
|
2759 for (int i = 0; i < nr; i++) |
|
2760 work[i] = 0.; |
|
2761 |
|
2762 for (int j = 0; j < nr; j++) |
|
2763 { |
|
2764 work[q_perm[j]] = 1.; |
|
2765 |
|
2766 for (int k = 0; k < nr; k++) |
|
2767 { |
|
2768 int iidx = q_perm[k]; |
|
2769 |
|
2770 if (work[iidx] != 0.) |
|
2771 { |
|
2772 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2773 work[iidx] = tmp; |
|
2774 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2775 { |
|
2776 int idx2 = q_perm[ridx(i)]; |
|
2777 work[idx2] = work[idx2] - tmp * data(i); |
|
2778 } |
|
2779 } |
|
2780 } |
|
2781 double atmp = 0; |
|
2782 for (int i = 0; i < j+1; i++) |
|
2783 { |
|
2784 atmp += ::abs(work[i]); |
|
2785 work[i] = 0.; |
|
2786 } |
|
2787 if (atmp > ainvnorm) |
|
2788 ainvnorm = atmp; |
|
2789 } |
|
2790 } |
|
2791 else |
|
2792 { |
|
2793 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2794 |
|
2795 for (int j = 0; j < b_nc; j++) |
|
2796 { |
|
2797 for (int i = 0; i < nr; i++) |
|
2798 work[i] = 0.; |
|
2799 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2800 work[b.ridx(i)] = b.data(i); |
|
2801 |
|
2802 for (int k = 0; k < nr; k++) |
|
2803 { |
|
2804 if (work[k] != 0.) |
|
2805 { |
|
2806 if (ridx(cidx(k)) != k) |
|
2807 { |
|
2808 err = -2; |
|
2809 goto triangular_error; |
|
2810 } |
|
2811 |
|
2812 Complex tmp = work[k] / data(cidx(k)); |
|
2813 work[k] = tmp; |
|
2814 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2815 { |
|
2816 int iidx = ridx(i); |
|
2817 work[iidx] = work[iidx] - tmp * data(i); |
|
2818 } |
|
2819 } |
|
2820 } |
|
2821 |
|
2822 // Count non-zeros in work vector and adjust space in |
|
2823 // retval if needed |
|
2824 int new_nnz = 0; |
|
2825 for (int i = 0; i < nr; i++) |
|
2826 if (work[i] != 0.) |
|
2827 new_nnz++; |
|
2828 |
|
2829 if (ii + new_nnz > x_nz) |
|
2830 { |
|
2831 // Resize the sparse matrix |
|
2832 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2833 retval.change_capacity (sz); |
|
2834 x_nz = sz; |
|
2835 } |
|
2836 |
|
2837 for (int i = 0; i < nr; i++) |
|
2838 if (work[i] != 0.) |
|
2839 { |
|
2840 retval.xridx(ii) = i; |
|
2841 retval.xdata(ii++) = work[i]; |
|
2842 } |
|
2843 retval.xcidx(j+1) = ii; |
|
2844 } |
|
2845 |
|
2846 retval.maybe_compress (); |
|
2847 |
|
2848 // Calculation of 1-norm of inv(*this) |
|
2849 for (int i = 0; i < nr; i++) |
|
2850 work[i] = 0.; |
|
2851 |
|
2852 for (int j = 0; j < nr; j++) |
|
2853 { |
|
2854 work[j] = 1.; |
|
2855 |
|
2856 for (int k = j; k < nr; k++) |
|
2857 { |
|
2858 |
|
2859 if (work[k] != 0.) |
|
2860 { |
|
2861 Complex tmp = work[k] / data(cidx(k)); |
|
2862 work[k] = tmp; |
|
2863 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2864 { |
|
2865 int iidx = ridx(i); |
|
2866 work[iidx] = work[iidx] - tmp * data(i); |
|
2867 } |
|
2868 } |
|
2869 } |
|
2870 double atmp = 0; |
|
2871 for (int i = j; i < nr; i++) |
|
2872 { |
|
2873 atmp += ::abs(work[i]); |
|
2874 work[i] = 0.; |
|
2875 } |
|
2876 if (atmp > ainvnorm) |
|
2877 ainvnorm = atmp; |
|
2878 } |
|
2879 |
|
2880 } |
|
2881 |
|
2882 rcond = 1. / ainvnorm / anorm; |
|
2883 |
|
2884 triangular_error: |
|
2885 if (err != 0) |
|
2886 { |
|
2887 if (sing_handler) |
|
2888 sing_handler (rcond); |
|
2889 else |
|
2890 (*current_liboctave_error_handler) |
|
2891 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2892 rcond); |
|
2893 } |
|
2894 |
|
2895 volatile double rcond_plus_one = rcond + 1.0; |
|
2896 |
|
2897 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2898 { |
|
2899 err = -2; |
|
2900 |
|
2901 if (sing_handler) |
|
2902 sing_handler (rcond); |
|
2903 else |
|
2904 (*current_liboctave_error_handler) |
|
2905 ("matrix singular to machine precision, rcond = %g", |
|
2906 rcond); |
|
2907 } |
|
2908 } |
|
2909 else |
|
2910 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2911 } |
|
2912 |
|
2913 return retval; |
|
2914 } |
|
2915 |
|
2916 ComplexMatrix |
|
2917 SparseComplexMatrix::trisolve (SparseType &mattype, const Matrix& b, int& err, |
|
2918 double& rcond, |
|
2919 solve_singularity_handler sing_handler) const |
|
2920 { |
|
2921 ComplexMatrix retval; |
|
2922 |
|
2923 int nr = rows (); |
|
2924 int nc = cols (); |
|
2925 err = 0; |
|
2926 |
|
2927 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2928 (*current_liboctave_error_handler) |
|
2929 ("matrix dimension mismatch solution of linear equations"); |
|
2930 else |
|
2931 { |
|
2932 // Print spparms("spumoni") info if requested |
|
2933 volatile int typ = mattype.type (); |
|
2934 mattype.info (); |
|
2935 |
|
2936 if (typ == SparseType::Tridiagonal_Hermitian) |
|
2937 { |
|
2938 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
2939 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
2940 |
|
2941 if (mattype.is_dense ()) |
|
2942 { |
|
2943 int ii = 0; |
|
2944 |
|
2945 for (int j = 0; j < nc-1; j++) |
|
2946 { |
|
2947 D[j] = data(ii++); |
|
2948 DL[j] = data(ii); |
|
2949 ii += 2; |
|
2950 } |
|
2951 D[nc-1] = data(ii); |
|
2952 } |
|
2953 else |
|
2954 { |
|
2955 D[0] = 0.; |
|
2956 for (int i = 0; i < nr - 1; i++) |
|
2957 { |
|
2958 D[i+1] = 0.; |
|
2959 DL[i] = 0.; |
|
2960 } |
|
2961 |
|
2962 for (int j = 0; j < nc; j++) |
|
2963 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2964 { |
|
2965 if (ridx(i) == j) |
|
2966 D[j] = data(i); |
|
2967 else if (ridx(i) == j + 1) |
|
2968 DL[j] = data(i); |
|
2969 } |
|
2970 } |
|
2971 |
|
2972 int b_nc = b.cols(); |
|
2973 retval = ComplexMatrix (b); |
|
2974 Complex *result = retval.fortran_vec (); |
|
2975 |
|
2976 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, |
|
2977 b.rows(), err)); |
|
2978 |
|
2979 if (f77_exception_encountered) |
|
2980 (*current_liboctave_error_handler) |
|
2981 ("unrecoverable error in zptsv"); |
|
2982 else if (err != 0) |
|
2983 { |
|
2984 err = 0; |
|
2985 mattype.mark_as_unsymmetric (); |
|
2986 typ = SparseType::Tridiagonal; |
|
2987 } |
|
2988 else |
|
2989 rcond = 1.; |
|
2990 } |
|
2991 |
|
2992 if (typ == SparseType::Tridiagonal) |
|
2993 { |
|
2994 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); |
|
2995 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
2996 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
2997 |
|
2998 if (mattype.is_dense ()) |
|
2999 { |
|
3000 int ii = 0; |
|
3001 |
|
3002 for (int j = 0; j < nc-1; j++) |
|
3003 { |
|
3004 D[j] = data(ii++); |
|
3005 DL[j] = data(ii++); |
|
3006 DU[j] = data(ii++); |
|
3007 } |
|
3008 D[nc-1] = data(ii); |
|
3009 } |
|
3010 else |
|
3011 { |
|
3012 D[0] = 0.; |
|
3013 for (int i = 0; i < nr - 1; i++) |
|
3014 { |
|
3015 D[i+1] = 0.; |
|
3016 DL[i] = 0.; |
|
3017 DU[i] = 0.; |
|
3018 } |
|
3019 |
|
3020 for (int j = 0; j < nc; j++) |
|
3021 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3022 { |
|
3023 if (ridx(i) == j) |
|
3024 D[j] = data(i); |
|
3025 else if (ridx(i) == j + 1) |
|
3026 DL[j] = data(i); |
|
3027 else if (ridx(i) == j - 1) |
|
3028 DU[j] = data(i); |
|
3029 } |
|
3030 } |
|
3031 |
|
3032 int b_nc = b.cols(); |
|
3033 retval = ComplexMatrix (b); |
|
3034 Complex *result = retval.fortran_vec (); |
|
3035 |
|
3036 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, |
|
3037 b.rows(), err)); |
|
3038 |
|
3039 if (f77_exception_encountered) |
|
3040 (*current_liboctave_error_handler) |
|
3041 ("unrecoverable error in zgtsv"); |
|
3042 else if (err != 0) |
|
3043 { |
|
3044 rcond = 0.; |
|
3045 err = -2; |
|
3046 |
|
3047 if (sing_handler) |
|
3048 sing_handler (rcond); |
|
3049 else |
|
3050 (*current_liboctave_error_handler) |
|
3051 ("matrix singular to machine precision"); |
|
3052 |
|
3053 } |
|
3054 else |
|
3055 rcond = 1.; |
|
3056 } |
|
3057 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3058 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3059 } |
|
3060 |
|
3061 return retval; |
|
3062 } |
|
3063 |
|
3064 SparseComplexMatrix |
|
3065 SparseComplexMatrix::trisolve (SparseType &mattype, const SparseMatrix& b, |
|
3066 int& err, double& rcond, |
|
3067 solve_singularity_handler sing_handler) const |
|
3068 { |
|
3069 SparseComplexMatrix retval; |
|
3070 |
|
3071 int nr = rows (); |
|
3072 int nc = cols (); |
|
3073 err = 0; |
|
3074 |
|
3075 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3076 (*current_liboctave_error_handler) |
|
3077 ("matrix dimension mismatch solution of linear equations"); |
|
3078 else |
|
3079 { |
|
3080 // Print spparms("spumoni") info if requested |
|
3081 int typ = mattype.type (); |
|
3082 mattype.info (); |
|
3083 |
|
3084 // Note can't treat symmetric case as there is no dpttrf function |
|
3085 if (typ == SparseType::Tridiagonal || |
|
3086 typ == SparseType::Tridiagonal_Hermitian) |
|
3087 { |
|
3088 OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2); |
|
3089 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); |
|
3090 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
3091 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
3092 Array<int> ipvt (nr); |
|
3093 int *pipvt = ipvt.fortran_vec (); |
|
3094 |
|
3095 if (mattype.is_dense ()) |
|
3096 { |
|
3097 int ii = 0; |
|
3098 |
|
3099 for (int j = 0; j < nc-1; j++) |
|
3100 { |
|
3101 D[j] = data(ii++); |
|
3102 DL[j] = data(ii++); |
|
3103 DU[j] = data(ii++); |
|
3104 } |
|
3105 D[nc-1] = data(ii); |
|
3106 } |
|
3107 else |
|
3108 { |
|
3109 D[0] = 0.; |
|
3110 for (int i = 0; i < nr - 1; i++) |
|
3111 { |
|
3112 D[i+1] = 0.; |
|
3113 DL[i] = 0.; |
|
3114 DU[i] = 0.; |
|
3115 } |
|
3116 |
|
3117 for (int j = 0; j < nc; j++) |
|
3118 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3119 { |
|
3120 if (ridx(i) == j) |
|
3121 D[j] = data(i); |
|
3122 else if (ridx(i) == j + 1) |
|
3123 DL[j] = data(i); |
|
3124 else if (ridx(i) == j - 1) |
|
3125 DU[j] = data(i); |
|
3126 } |
|
3127 } |
|
3128 |
|
3129 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); |
|
3130 |
|
3131 if (f77_exception_encountered) |
|
3132 (*current_liboctave_error_handler) |
|
3133 ("unrecoverable error in zgttrf"); |
|
3134 else |
|
3135 { |
|
3136 rcond = 0.0; |
|
3137 if (err != 0) |
|
3138 { |
|
3139 err = -2; |
|
3140 |
|
3141 if (sing_handler) |
|
3142 sing_handler (rcond); |
|
3143 else |
|
3144 (*current_liboctave_error_handler) |
|
3145 ("matrix singular to machine precision"); |
|
3146 |
|
3147 } |
|
3148 else |
|
3149 { |
|
3150 char job = 'N'; |
|
3151 volatile int x_nz = b.nnz (); |
|
3152 int b_nc = b.cols (); |
|
3153 retval = SparseComplexMatrix (nr, b_nc, x_nz); |
|
3154 retval.xcidx(0) = 0; |
|
3155 volatile int ii = 0; |
|
3156 |
|
3157 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
3158 |
|
3159 for (volatile int j = 0; j < b_nc; j++) |
|
3160 { |
|
3161 for (int i = 0; i < nr; i++) |
|
3162 work[i] = 0.; |
|
3163 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
3164 work[b.ridx(i)] = b.data(i); |
|
3165 |
|
3166 F77_XFCN (zgttrs, ZGTTRS, |
|
3167 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3168 nr, 1, DL, D, DU, DU2, pipvt, |
|
3169 work, b.rows (), err |
|
3170 F77_CHAR_ARG_LEN (1))); |
|
3171 |
|
3172 if (f77_exception_encountered) |
|
3173 { |
|
3174 (*current_liboctave_error_handler) |
|
3175 ("unrecoverable error in zgttrs"); |
|
3176 break; |
|
3177 } |
|
3178 |
|
3179 // Count non-zeros in work vector and adjust |
|
3180 // space in retval if needed |
|
3181 int new_nnz = 0; |
|
3182 for (int i = 0; i < nr; i++) |
|
3183 if (work[i] != 0.) |
|
3184 new_nnz++; |
|
3185 |
|
3186 if (ii + new_nnz > x_nz) |
|
3187 { |
|
3188 // Resize the sparse matrix |
|
3189 int sz = new_nnz * (b_nc - j) + x_nz; |
|
3190 retval.change_capacity (sz); |
|
3191 x_nz = sz; |
|
3192 } |
|
3193 |
|
3194 for (int i = 0; i < nr; i++) |
|
3195 if (work[i] != 0.) |
|
3196 { |
|
3197 retval.xridx(ii) = i; |
|
3198 retval.xdata(ii++) = work[i]; |
|
3199 } |
|
3200 retval.xcidx(j+1) = ii; |
|
3201 } |
|
3202 |
|
3203 retval.maybe_compress (); |
|
3204 } |
|
3205 } |
|
3206 } |
|
3207 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3208 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3209 } |
|
3210 |
|
3211 return retval; |
|
3212 } |
|
3213 |
|
3214 ComplexMatrix |
|
3215 SparseComplexMatrix::trisolve (SparseType &mattype, const ComplexMatrix& b, |
|
3216 int& err, double& rcond, |
|
3217 solve_singularity_handler sing_handler) const |
|
3218 { |
|
3219 ComplexMatrix retval; |
|
3220 |
|
3221 int nr = rows (); |
|
3222 int nc = cols (); |
|
3223 err = 0; |
|
3224 |
|
3225 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3226 (*current_liboctave_error_handler) |
|
3227 ("matrix dimension mismatch solution of linear equations"); |
|
3228 else |
|
3229 { |
|
3230 // Print spparms("spumoni") info if requested |
|
3231 volatile int typ = mattype.type (); |
|
3232 mattype.info (); |
|
3233 |
|
3234 // Note can't treat symmetric case as there is no dpttrf function |
|
3235 if (typ == SparseType::Tridiagonal_Hermitian) |
|
3236 { |
|
3237 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
3238 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
3239 |
|
3240 if (mattype.is_dense ()) |
|
3241 { |
|
3242 int ii = 0; |
|
3243 |
|
3244 for (int j = 0; j < nc-1; j++) |
|
3245 { |
|
3246 D[j] = data(ii++); |
|
3247 DL[j] = data(ii); |
|
3248 ii += 2; |
|
3249 } |
|
3250 D[nc-1] = data(ii); |
|
3251 } |
|
3252 else |
|
3253 { |
|
3254 D[0] = 0.; |
|
3255 for (int i = 0; i < nr - 1; i++) |
|
3256 { |
|
3257 D[i+1] = 0.; |
|
3258 DL[i] = 0.; |
|
3259 } |
|
3260 |
|
3261 for (int j = 0; j < nc; j++) |
|
3262 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3263 { |
|
3264 if (ridx(i) == j) |
|
3265 D[j] = data(i); |
|
3266 else if (ridx(i) == j + 1) |
|
3267 DL[j] = data(i); |
|
3268 } |
|
3269 } |
|
3270 |
|
3271 int b_nr = b.rows (); |
|
3272 int b_nc = b.cols(); |
|
3273 rcond = 1.; |
|
3274 |
|
3275 retval = ComplexMatrix (b); |
|
3276 Complex *result = retval.fortran_vec (); |
|
3277 |
|
3278 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, |
|
3279 b_nr, err)); |
|
3280 |
|
3281 if (f77_exception_encountered) |
|
3282 { |
|
3283 (*current_liboctave_error_handler) |
|
3284 ("unrecoverable error in zptsv"); |
|
3285 err = -1; |
|
3286 } |
|
3287 else if (err != 0) |
|
3288 { |
|
3289 err = 0; |
|
3290 mattype.mark_as_unsymmetric (); |
|
3291 typ = SparseType::Tridiagonal; |
|
3292 } |
|
3293 } |
|
3294 |
|
3295 if (typ == SparseType::Tridiagonal) |
|
3296 { |
|
3297 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); |
|
3298 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
3299 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
3300 |
|
3301 if (mattype.is_dense ()) |
|
3302 { |
|
3303 int ii = 0; |
|
3304 |
|
3305 for (int j = 0; j < nc-1; j++) |
|
3306 { |
|
3307 D[j] = data(ii++); |
|
3308 DL[j] = data(ii++); |
|
3309 DU[j] = data(ii++); |
|
3310 } |
|
3311 D[nc-1] = data(ii); |
|
3312 } |
|
3313 else |
|
3314 { |
|
3315 D[0] = 0.; |
|
3316 for (int i = 0; i < nr - 1; i++) |
|
3317 { |
|
3318 D[i+1] = 0.; |
|
3319 DL[i] = 0.; |
|
3320 DU[i] = 0.; |
|
3321 } |
|
3322 |
|
3323 for (int j = 0; j < nc; j++) |
|
3324 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3325 { |
|
3326 if (ridx(i) == j) |
|
3327 D[j] = data(i); |
|
3328 else if (ridx(i) == j + 1) |
|
3329 DL[j] = data(i); |
|
3330 else if (ridx(i) == j - 1) |
|
3331 DU[j] = data(i); |
|
3332 } |
|
3333 } |
|
3334 |
|
3335 int b_nr = b.rows(); |
|
3336 int b_nc = b.cols(); |
|
3337 rcond = 1.; |
|
3338 |
|
3339 retval = ComplexMatrix (b); |
|
3340 Complex *result = retval.fortran_vec (); |
|
3341 |
|
3342 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, |
|
3343 b_nr, err)); |
|
3344 |
|
3345 if (f77_exception_encountered) |
|
3346 { |
|
3347 (*current_liboctave_error_handler) |
|
3348 ("unrecoverable error in zgtsv"); |
|
3349 err = -1; |
|
3350 } |
|
3351 else if (err != 0) |
|
3352 { |
|
3353 rcond = 0.; |
|
3354 err = -2; |
|
3355 |
|
3356 if (sing_handler) |
|
3357 sing_handler (rcond); |
|
3358 else |
|
3359 (*current_liboctave_error_handler) |
|
3360 ("matrix singular to machine precision"); |
|
3361 } |
|
3362 } |
|
3363 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3364 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3365 } |
|
3366 |
|
3367 return retval; |
|
3368 } |
|
3369 |
|
3370 SparseComplexMatrix |
|
3371 SparseComplexMatrix::trisolve (SparseType &mattype, |
|
3372 const SparseComplexMatrix& b, int& err, double& rcond, |
|
3373 solve_singularity_handler sing_handler) const |
|
3374 { |
|
3375 SparseComplexMatrix retval; |
|
3376 |
|
3377 int nr = rows (); |
|
3378 int nc = cols (); |
|
3379 err = 0; |
|
3380 |
|
3381 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3382 (*current_liboctave_error_handler) |
|
3383 ("matrix dimension mismatch solution of linear equations"); |
|
3384 else |
|
3385 { |
|
3386 // Print spparms("spumoni") info if requested |
|
3387 int typ = mattype.type (); |
|
3388 mattype.info (); |
|
3389 |
|
3390 // Note can't treat symmetric case as there is no dpttrf function |
|
3391 if (typ == SparseType::Tridiagonal || |
|
3392 typ == SparseType::Tridiagonal_Hermitian) |
|
3393 { |
|
3394 OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2); |
|
3395 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); |
|
3396 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
3397 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
3398 Array<int> ipvt (nr); |
|
3399 int *pipvt = ipvt.fortran_vec (); |
|
3400 |
|
3401 if (mattype.is_dense ()) |
|
3402 { |
|
3403 int ii = 0; |
|
3404 |
|
3405 for (int j = 0; j < nc-1; j++) |
|
3406 { |
|
3407 D[j] = data(ii++); |
|
3408 DL[j] = data(ii++); |
|
3409 DU[j] = data(ii++); |
|
3410 } |
|
3411 D[nc-1] = data(ii); |
|
3412 } |
|
3413 else |
|
3414 { |
|
3415 D[0] = 0.; |
|
3416 for (int i = 0; i < nr - 1; i++) |
|
3417 { |
|
3418 D[i+1] = 0.; |
|
3419 DL[i] = 0.; |
|
3420 DU[i] = 0.; |
|
3421 } |
|
3422 |
|
3423 for (int j = 0; j < nc; j++) |
|
3424 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3425 { |
|
3426 if (ridx(i) == j) |
|
3427 D[j] = data(i); |
|
3428 else if (ridx(i) == j + 1) |
|
3429 DL[j] = data(i); |
|
3430 else if (ridx(i) == j - 1) |
|
3431 DU[j] = data(i); |
|
3432 } |
|
3433 } |
|
3434 |
|
3435 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); |
|
3436 |
|
3437 if (f77_exception_encountered) |
|
3438 (*current_liboctave_error_handler) |
|
3439 ("unrecoverable error in zgttrf"); |
|
3440 else |
|
3441 { |
|
3442 rcond = 0.0; |
|
3443 if (err != 0) |
|
3444 { |
|
3445 err = -2; |
|
3446 |
|
3447 if (sing_handler) |
|
3448 sing_handler (rcond); |
|
3449 else |
|
3450 (*current_liboctave_error_handler) |
|
3451 ("matrix singular to machine precision"); |
|
3452 } |
|
3453 else |
|
3454 { |
|
3455 rcond = 1.; |
|
3456 char job = 'N'; |
|
3457 int b_nr = b.rows (); |
|
3458 int b_nc = b.cols (); |
|
3459 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); |
|
3460 |
|
3461 // Take a first guess that the number of non-zero terms |
|
3462 // will be as many as in b |
|
3463 volatile int x_nz = b.nnz (); |
|
3464 volatile int ii = 0; |
|
3465 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
3466 |
|
3467 retval.xcidx(0) = 0; |
|
3468 for (volatile int j = 0; j < b_nc; j++) |
|
3469 { |
|
3470 |
|
3471 for (int i = 0; i < b_nr; i++) |
|
3472 Bx[i] = b (i,j); |
|
3473 |
|
3474 F77_XFCN (zgttrs, ZGTTRS, |
|
3475 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3476 nr, 1, DL, D, DU, DU2, pipvt, |
|
3477 Bx, b_nr, err |
|
3478 F77_CHAR_ARG_LEN (1))); |
|
3479 |
|
3480 if (f77_exception_encountered) |
|
3481 { |
|
3482 (*current_liboctave_error_handler) |
|
3483 ("unrecoverable error in zgttrs"); |
|
3484 break; |
|
3485 } |
|
3486 |
|
3487 if (err != 0) |
|
3488 { |
|
3489 (*current_liboctave_error_handler) |
|
3490 ("SparseComplexMatrix::solve solve failed"); |
|
3491 |
|
3492 err = -1; |
|
3493 break; |
|
3494 } |
|
3495 |
|
3496 // Count non-zeros in work vector and adjust |
|
3497 // space in retval if needed |
|
3498 int new_nnz = 0; |
|
3499 for (int i = 0; i < nr; i++) |
|
3500 if (Bx[i] != 0.) |
|
3501 new_nnz++; |
|
3502 |
|
3503 if (ii + new_nnz > x_nz) |
|
3504 { |
|
3505 // Resize the sparse matrix |
|
3506 int sz = new_nnz * (b_nc - j) + x_nz; |
|
3507 retval.change_capacity (sz); |
|
3508 x_nz = sz; |
|
3509 } |
|
3510 |
|
3511 for (int i = 0; i < nr; i++) |
|
3512 if (Bx[i] != 0.) |
|
3513 { |
|
3514 retval.xridx(ii) = i; |
|
3515 retval.xdata(ii++) = Bx[i]; |
|
3516 } |
|
3517 |
|
3518 retval.xcidx(j+1) = ii; |
|
3519 } |
|
3520 |
|
3521 retval.maybe_compress (); |
|
3522 } |
|
3523 } |
|
3524 } |
|
3525 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3526 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3527 } |
|
3528 |
|
3529 return retval; |
|
3530 } |
|
3531 |
|
3532 ComplexMatrix |
|
3533 SparseComplexMatrix::bsolve (SparseType &mattype, const Matrix& b, int& err, |
|
3534 double& rcond, |
|
3535 solve_singularity_handler sing_handler) const |
|
3536 { |
|
3537 ComplexMatrix retval; |
|
3538 |
|
3539 int nr = rows (); |
|
3540 int nc = cols (); |
|
3541 err = 0; |
|
3542 |
|
3543 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3544 (*current_liboctave_error_handler) |
|
3545 ("matrix dimension mismatch solution of linear equations"); |
|
3546 else |
|
3547 { |
|
3548 // Print spparms("spumoni") info if requested |
|
3549 volatile int typ = mattype.type (); |
|
3550 mattype.info (); |
|
3551 |
|
3552 if (typ == SparseType::Banded_Hermitian) |
|
3553 { |
|
3554 int n_lower = mattype.nlower (); |
|
3555 int ldm = n_lower + 1; |
|
3556 ComplexMatrix m_band (ldm, nc); |
|
3557 Complex *tmp_data = m_band.fortran_vec (); |
|
3558 |
|
3559 if (! mattype.is_dense ()) |
|
3560 { |
|
3561 int ii = 0; |
|
3562 |
|
3563 for (int j = 0; j < ldm; j++) |
|
3564 for (int i = 0; i < nc; i++) |
|
3565 tmp_data[ii++] = 0.; |
|
3566 } |
|
3567 |
|
3568 for (int j = 0; j < nc; j++) |
|
3569 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3570 { |
|
3571 int ri = ridx (i); |
|
3572 if (ri >= j) |
|
3573 m_band(ri - j, j) = data(i); |
|
3574 } |
|
3575 |
|
3576 // Calculate the norm of the matrix, for later use. |
|
3577 // double anorm = m_band.abs().sum().row(0).max(); |
|
3578 |
|
3579 char job = 'L'; |
|
3580 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3581 nr, n_lower, tmp_data, ldm, err |
|
3582 F77_CHAR_ARG_LEN (1))); |
|
3583 |
|
3584 if (f77_exception_encountered) |
|
3585 (*current_liboctave_error_handler) |
|
3586 ("unrecoverable error in zpbtrf"); |
|
3587 else |
|
3588 { |
|
3589 rcond = 0.0; |
|
3590 if (err != 0) |
|
3591 { |
|
3592 // Matrix is not positive definite!! Fall through to |
|
3593 // unsymmetric banded solver. |
|
3594 mattype.mark_as_unsymmetric (); |
|
3595 typ = SparseType::Banded; |
|
3596 err = 0; |
|
3597 } |
|
3598 else |
|
3599 { |
|
3600 // Unfortunately, the time to calculate the condition |
|
3601 // number is dominant for narrow banded matrices and |
|
3602 // so we rely on the "err" flag from xPBTRF to flag |
|
3603 // singularity. The commented code below is left here |
|
3604 // for reference |
|
3605 |
|
3606 //Array<double> z (3 * nr); |
|
3607 //Complex *pz = z.fortran_vec (); |
|
3608 //Array<int> iz (nr); |
|
3609 //int *piz = iz.fortran_vec (); |
|
3610 // |
|
3611 //F77_XFCN (zpbcon, ZGBCON, |
|
3612 // (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3613 // nr, n_lower, tmp_data, ldm, |
|
3614 // anorm, rcond, pz, piz, err |
|
3615 // F77_CHAR_ARG_LEN (1))); |
|
3616 // |
|
3617 // |
|
3618 //if (f77_exception_encountered) |
|
3619 // (*current_liboctave_error_handler) |
|
3620 // ("unrecoverable error in zpbcon"); |
|
3621 // |
|
3622 //if (err != 0) |
|
3623 // err = -2; |
|
3624 // |
|
3625 //volatile double rcond_plus_one = rcond + 1.0; |
|
3626 // |
|
3627 //if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
3628 // { |
|
3629 // err = -2; |
|
3630 // |
|
3631 // if (sing_handler) |
|
3632 // sing_handler (rcond); |
|
3633 // else |
|
3634 // (*current_liboctave_error_handler) |
|
3635 // ("matrix singular to machine precision, rcond = %g", |
|
3636 // rcond); |
|
3637 // } |
|
3638 //else |
|
3639 // REST OF CODE, EXCEPT rcond=1 |
|
3640 |
|
3641 rcond = 1.; |
|
3642 retval = ComplexMatrix (b); |
|
3643 Complex *result = retval.fortran_vec (); |
|
3644 |
|
3645 int b_nc = b.cols (); |
|
3646 |
|
3647 F77_XFCN (zpbtrs, ZPBTRS, |
|
3648 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3649 nr, n_lower, b_nc, tmp_data, |
|
3650 ldm, result, b.rows(), err |
|
3651 F77_CHAR_ARG_LEN (1))); |
|
3652 |
|
3653 if (f77_exception_encountered) |
|
3654 (*current_liboctave_error_handler) |
|
3655 ("unrecoverable error in zpbtrs"); |
|
3656 |
|
3657 if (err != 0) |
|
3658 { |
|
3659 (*current_liboctave_error_handler) |
|
3660 ("SparseMatrix::solve solve failed"); |
|
3661 err = -1; |
|
3662 } |
|
3663 } |
|
3664 } |
|
3665 } |
|
3666 |
|
3667 if (typ == SparseType::Banded) |
|
3668 { |
|
3669 // Create the storage for the banded form of the sparse matrix |
|
3670 int n_upper = mattype.nupper (); |
|
3671 int n_lower = mattype.nlower (); |
|
3672 int ldm = n_upper + 2 * n_lower + 1; |
|
3673 |
|
3674 ComplexMatrix m_band (ldm, nc); |
|
3675 Complex *tmp_data = m_band.fortran_vec (); |
|
3676 |
|
3677 if (! mattype.is_dense ()) |
|
3678 { |
|
3679 int ii = 0; |
|
3680 |
|
3681 for (int j = 0; j < ldm; j++) |
|
3682 for (int i = 0; i < nc; i++) |
|
3683 tmp_data[ii++] = 0.; |
|
3684 } |
|
3685 |
|
3686 for (int j = 0; j < nc; j++) |
|
3687 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3688 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
3689 |
|
3690 Array<int> ipvt (nr); |
|
3691 int *pipvt = ipvt.fortran_vec (); |
|
3692 |
|
3693 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
3694 ldm, pipvt, err)); |
|
3695 |
|
3696 if (f77_exception_encountered) |
|
3697 (*current_liboctave_error_handler) |
|
3698 ("unrecoverable error in zgbtrf"); |
|
3699 else |
|
3700 { |
|
3701 // Throw-away extra info LAPACK gives so as to not |
|
3702 // change output. |
|
3703 rcond = 0.0; |
|
3704 if (err != 0) |
|
3705 { |
|
3706 err = -2; |
|
3707 |
|
3708 if (sing_handler) |
|
3709 sing_handler (rcond); |
|
3710 else |
|
3711 (*current_liboctave_error_handler) |
|
3712 ("matrix singular to machine precision"); |
|
3713 |
|
3714 } |
|
3715 else |
|
3716 { |
|
3717 char job = '1'; |
|
3718 |
|
3719 // Unfortunately, the time to calculate the condition |
|
3720 // number is dominant for narrow banded matrices and |
|
3721 // so we rely on the "err" flag from xPBTRF to flag |
|
3722 // singularity. The commented code below is left here |
|
3723 // for reference |
|
3724 |
|
3725 //F77_XFCN (zgbcon, ZGBCON, |
|
3726 // (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3727 // nc, n_lower, n_upper, tmp_data, ldm, pipvt, |
|
3728 // anorm, rcond, pz, piz, err |
|
3729 // F77_CHAR_ARG_LEN (1))); |
|
3730 // |
|
3731 //if (f77_exception_encountered) |
|
3732 // (*current_liboctave_error_handler) |
|
3733 // ("unrecoverable error in zgbcon"); |
|
3734 // |
|
3735 // if (err != 0) |
|
3736 // err = -2; |
|
3737 // |
|
3738 //volatile double rcond_plus_one = rcond + 1.0; |
|
3739 // |
|
3740 //if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
3741 // { |
|
3742 // err = -2; |
|
3743 // |
|
3744 // if (sing_handler) |
|
3745 // sing_handler (rcond); |
|
3746 // else |
|
3747 // (*current_liboctave_error_handler) |
|
3748 // ("matrix singular to machine precision, rcond = %g", |
|
3749 // rcond); |
|
3750 // } |
|
3751 //else |
|
3752 // REST OF CODE, EXCEPT rcond=1 |
|
3753 |
|
3754 rcond = 1.; |
|
3755 retval = ComplexMatrix (b); |
|
3756 Complex *result = retval.fortran_vec (); |
|
3757 |
|
3758 int b_nc = b.cols (); |
|
3759 |
|
3760 job = 'N'; |
|
3761 F77_XFCN (zgbtrs, ZGBTRS, |
|
3762 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3763 nr, n_lower, n_upper, b_nc, tmp_data, |
|
3764 ldm, pipvt, result, b.rows(), err |
|
3765 F77_CHAR_ARG_LEN (1))); |
|
3766 |
|
3767 if (f77_exception_encountered) |
|
3768 (*current_liboctave_error_handler) |
|
3769 ("unrecoverable error in zgbtrs"); |
|
3770 } |
|
3771 } |
|
3772 } |
|
3773 else if (typ != SparseType::Banded_Hermitian) |
|
3774 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3775 } |
|
3776 |
|
3777 return retval; |
|
3778 } |
|
3779 |
|
3780 SparseComplexMatrix |
|
3781 SparseComplexMatrix::bsolve (SparseType &mattype, const SparseMatrix& b, |
|
3782 int& err, double& rcond, |
|
3783 solve_singularity_handler sing_handler) const |
|
3784 { |
|
3785 SparseComplexMatrix retval; |
|
3786 |
|
3787 int nr = rows (); |
|
3788 int nc = cols (); |
|
3789 err = 0; |
|
3790 |
|
3791 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3792 (*current_liboctave_error_handler) |
|
3793 ("matrix dimension mismatch solution of linear equations"); |
|
3794 else |
|
3795 { |
|
3796 // Print spparms("spumoni") info if requested |
|
3797 volatile int typ = mattype.type (); |
|
3798 mattype.info (); |
|
3799 |
|
3800 if (typ == SparseType::Banded_Hermitian) |
|
3801 { |
|
3802 int n_lower = mattype.nlower (); |
|
3803 int ldm = n_lower + 1; |
|
3804 |
|
3805 ComplexMatrix m_band (ldm, nc); |
|
3806 Complex *tmp_data = m_band.fortran_vec (); |
|
3807 |
|
3808 if (! mattype.is_dense ()) |
|
3809 { |
|
3810 int ii = 0; |
|
3811 |
|
3812 for (int j = 0; j < ldm; j++) |
|
3813 for (int i = 0; i < nc; i++) |
|
3814 tmp_data[ii++] = 0.; |
|
3815 } |
|
3816 |
|
3817 for (int j = 0; j < nc; j++) |
|
3818 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3819 { |
|
3820 int ri = ridx (i); |
|
3821 if (ri >= j) |
|
3822 m_band(ri - j, j) = data(i); |
|
3823 } |
|
3824 |
|
3825 char job = 'L'; |
|
3826 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3827 nr, n_lower, tmp_data, ldm, err |
|
3828 F77_CHAR_ARG_LEN (1))); |
|
3829 |
|
3830 if (f77_exception_encountered) |
|
3831 (*current_liboctave_error_handler) |
|
3832 ("unrecoverable error in zpbtrf"); |
|
3833 else |
|
3834 { |
|
3835 rcond = 0.0; |
|
3836 if (err != 0) |
|
3837 { |
|
3838 mattype.mark_as_unsymmetric (); |
|
3839 typ = SparseType::Banded; |
|
3840 err = 0; |
|
3841 } |
|
3842 else |
|
3843 { |
|
3844 rcond = 1.; |
|
3845 int b_nr = b.rows (); |
|
3846 int b_nc = b.cols (); |
|
3847 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); |
|
3848 |
|
3849 // Take a first guess that the number of non-zero terms |
|
3850 // will be as many as in b |
|
3851 volatile int x_nz = b.nnz (); |
|
3852 volatile int ii = 0; |
|
3853 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
3854 |
|
3855 retval.xcidx(0) = 0; |
|
3856 for (volatile int j = 0; j < b_nc; j++) |
|
3857 { |
|
3858 for (int i = 0; i < b_nr; i++) |
|
3859 Bx[i] = b.elem (i, j); |
|
3860 |
|
3861 F77_XFCN (zpbtrs, ZPBTRS, |
|
3862 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3863 nr, n_lower, 1, tmp_data, |
|
3864 ldm, Bx, b_nr, err |
|
3865 F77_CHAR_ARG_LEN (1))); |
|
3866 |
|
3867 if (f77_exception_encountered) |
|
3868 { |
|
3869 (*current_liboctave_error_handler) |
|
3870 ("unrecoverable error in dpbtrs"); |
|
3871 err = -1; |
|
3872 break; |
|
3873 } |
|
3874 |
|
3875 if (err != 0) |
|
3876 { |
|
3877 (*current_liboctave_error_handler) |
|
3878 ("SparseComplexMatrix::solve solve failed"); |
|
3879 err = -1; |
|
3880 break; |
|
3881 } |
|
3882 |
|
3883 for (int i = 0; i < b_nr; i++) |
|
3884 { |
|
3885 Complex tmp = Bx[i]; |
|
3886 if (tmp != 0.0) |
|
3887 { |
|
3888 if (ii == x_nz) |
|
3889 { |
|
3890 // Resize the sparse matrix |
|
3891 int sz = x_nz * (b_nc - j) / b_nc; |
|
3892 sz = (sz > 10 ? sz : 10) + x_nz; |
|
3893 retval.change_capacity (sz); |
|
3894 x_nz = sz; |
|
3895 } |
|
3896 retval.xdata(ii) = tmp; |
|
3897 retval.xridx(ii++) = i; |
|
3898 } |
|
3899 } |
|
3900 retval.xcidx(j+1) = ii; |
|
3901 } |
|
3902 |
|
3903 retval.maybe_compress (); |
|
3904 } |
|
3905 } |
|
3906 } |
|
3907 |
|
3908 if (typ == SparseType::Banded) |
|
3909 { |
|
3910 // Create the storage for the banded form of the sparse matrix |
|
3911 int n_upper = mattype.nupper (); |
|
3912 int n_lower = mattype.nlower (); |
|
3913 int ldm = n_upper + 2 * n_lower + 1; |
|
3914 |
|
3915 ComplexMatrix m_band (ldm, nc); |
|
3916 Complex *tmp_data = m_band.fortran_vec (); |
|
3917 |
|
3918 if (! mattype.is_dense ()) |
|
3919 { |
|
3920 int ii = 0; |
|
3921 |
|
3922 for (int j = 0; j < ldm; j++) |
|
3923 for (int i = 0; i < nc; i++) |
|
3924 tmp_data[ii++] = 0.; |
|
3925 } |
|
3926 |
|
3927 for (int j = 0; j < nc; j++) |
|
3928 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3929 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
3930 |
|
3931 Array<int> ipvt (nr); |
|
3932 int *pipvt = ipvt.fortran_vec (); |
|
3933 |
|
3934 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
3935 ldm, pipvt, err)); |
|
3936 |
|
3937 if (f77_exception_encountered) |
|
3938 (*current_liboctave_error_handler) |
|
3939 ("unrecoverable error in zgbtrf"); |
|
3940 else |
|
3941 { |
|
3942 rcond = 0.0; |
|
3943 if (err != 0) |
|
3944 { |
|
3945 err = -2; |
|
3946 |
|
3947 if (sing_handler) |
|
3948 sing_handler (rcond); |
|
3949 else |
|
3950 (*current_liboctave_error_handler) |
|
3951 ("matrix singular to machine precision"); |
|
3952 |
|
3953 } |
|
3954 else |
|
3955 { |
|
3956 char job = 'N'; |
|
3957 volatile int x_nz = b.nnz (); |
|
3958 int b_nc = b.cols (); |
|
3959 retval = SparseComplexMatrix (nr, b_nc, x_nz); |
|
3960 retval.xcidx(0) = 0; |
|
3961 volatile int ii = 0; |
|
3962 |
|
3963 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
3964 |
|
3965 for (volatile int j = 0; j < b_nc; j++) |
|
3966 { |
|
3967 for (int i = 0; i < nr; i++) |
|
3968 work[i] = 0.; |
|
3969 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
3970 work[b.ridx(i)] = b.data(i); |
|
3971 |
|
3972 F77_XFCN (zgbtrs, ZGBTRS, |
|
3973 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3974 nr, n_lower, n_upper, 1, tmp_data, |
|
3975 ldm, pipvt, work, b.rows (), err |
|
3976 F77_CHAR_ARG_LEN (1))); |
|
3977 |
|
3978 if (f77_exception_encountered) |
|
3979 { |
|
3980 (*current_liboctave_error_handler) |
|
3981 ("unrecoverable error in zgbtrs"); |
|
3982 break; |
|
3983 } |
|
3984 |
|
3985 // Count non-zeros in work vector and adjust |
|
3986 // space in retval if needed |
|
3987 int new_nnz = 0; |
|
3988 for (int i = 0; i < nr; i++) |
|
3989 if (work[i] != 0.) |
|
3990 new_nnz++; |
|
3991 |
|
3992 if (ii + new_nnz > x_nz) |
|
3993 { |
|
3994 // Resize the sparse matrix |
|
3995 int sz = new_nnz * (b_nc - j) + x_nz; |
|
3996 retval.change_capacity (sz); |
|
3997 x_nz = sz; |
|
3998 } |
|
3999 |
|
4000 for (int i = 0; i < nr; i++) |
|
4001 if (work[i] != 0.) |
|
4002 { |
|
4003 retval.xridx(ii) = i; |
|
4004 retval.xdata(ii++) = work[i]; |
|
4005 } |
|
4006 retval.xcidx(j+1) = ii; |
|
4007 } |
|
4008 |
|
4009 retval.maybe_compress (); |
|
4010 } |
|
4011 } |
|
4012 } |
|
4013 else if (typ != SparseType::Banded_Hermitian) |
|
4014 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4015 } |
|
4016 |
|
4017 return retval; |
|
4018 } |
|
4019 |
|
4020 ComplexMatrix |
|
4021 SparseComplexMatrix::bsolve (SparseType &mattype, const ComplexMatrix& b, |
|
4022 int& err, double& rcond, |
|
4023 solve_singularity_handler sing_handler) const |
|
4024 { |
|
4025 ComplexMatrix retval; |
|
4026 |
|
4027 int nr = rows (); |
|
4028 int nc = cols (); |
|
4029 err = 0; |
|
4030 |
|
4031 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4032 (*current_liboctave_error_handler) |
|
4033 ("matrix dimension mismatch solution of linear equations"); |
|
4034 else |
|
4035 { |
|
4036 // Print spparms("spumoni") info if requested |
|
4037 volatile int typ = mattype.type (); |
|
4038 mattype.info (); |
|
4039 |
|
4040 if (typ == SparseType::Banded_Hermitian) |
|
4041 { |
|
4042 int n_lower = mattype.nlower (); |
|
4043 int ldm = n_lower + 1; |
|
4044 |
|
4045 ComplexMatrix m_band (ldm, nc); |
|
4046 Complex *tmp_data = m_band.fortran_vec (); |
|
4047 |
|
4048 if (! mattype.is_dense ()) |
|
4049 { |
|
4050 int ii = 0; |
|
4051 |
|
4052 for (int j = 0; j < ldm; j++) |
|
4053 for (int i = 0; i < nc; i++) |
|
4054 tmp_data[ii++] = 0.; |
|
4055 } |
|
4056 |
|
4057 for (int j = 0; j < nc; j++) |
|
4058 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4059 { |
|
4060 int ri = ridx (i); |
|
4061 if (ri >= j) |
|
4062 m_band(ri - j, j) = data(i); |
|
4063 } |
|
4064 |
|
4065 char job = 'L'; |
|
4066 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4067 nr, n_lower, tmp_data, ldm, err |
|
4068 F77_CHAR_ARG_LEN (1))); |
|
4069 |
|
4070 if (f77_exception_encountered) |
|
4071 (*current_liboctave_error_handler) |
|
4072 ("unrecoverable error in zpbtrf"); |
|
4073 else |
|
4074 { |
|
4075 rcond = 0.0; |
|
4076 if (err != 0) |
|
4077 { |
|
4078 // Matrix is not positive definite!! Fall through to |
|
4079 // unsymmetric banded solver. |
|
4080 mattype.mark_as_unsymmetric (); |
|
4081 typ = SparseType::Banded; |
|
4082 err = 0; |
|
4083 } |
|
4084 else |
|
4085 { |
|
4086 rcond = 1.; |
|
4087 int b_nr = b.rows (); |
|
4088 int b_nc = b.cols (); |
|
4089 retval = ComplexMatrix (b); |
|
4090 Complex *result = retval.fortran_vec (); |
|
4091 |
|
4092 F77_XFCN (zpbtrs, ZPBTRS, |
|
4093 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4094 nr, n_lower, b_nc, tmp_data, |
|
4095 ldm, result, b_nr, err |
|
4096 F77_CHAR_ARG_LEN (1))); |
|
4097 |
|
4098 if (f77_exception_encountered) |
|
4099 { |
|
4100 (*current_liboctave_error_handler) |
|
4101 ("unrecoverable error in zpbtrs"); |
|
4102 err = -1; |
|
4103 } |
|
4104 |
|
4105 if (err != 0) |
|
4106 { |
|
4107 (*current_liboctave_error_handler) |
|
4108 ("SparseComplexMatrix::solve solve failed"); |
|
4109 err = -1; |
|
4110 } |
|
4111 } |
|
4112 } |
|
4113 } |
|
4114 |
|
4115 if (typ == SparseType::Banded) |
|
4116 { |
|
4117 // Create the storage for the banded form of the sparse matrix |
|
4118 int n_upper = mattype.nupper (); |
|
4119 int n_lower = mattype.nlower (); |
|
4120 int ldm = n_upper + 2 * n_lower + 1; |
|
4121 |
|
4122 ComplexMatrix m_band (ldm, nc); |
|
4123 Complex *tmp_data = m_band.fortran_vec (); |
|
4124 |
|
4125 if (! mattype.is_dense ()) |
|
4126 { |
|
4127 int ii = 0; |
|
4128 |
|
4129 for (int j = 0; j < ldm; j++) |
|
4130 for (int i = 0; i < nc; i++) |
|
4131 tmp_data[ii++] = 0.; |
|
4132 } |
|
4133 |
|
4134 for (int j = 0; j < nc; j++) |
|
4135 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4136 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
4137 |
|
4138 Array<int> ipvt (nr); |
|
4139 int *pipvt = ipvt.fortran_vec (); |
|
4140 |
|
4141 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
4142 ldm, pipvt, err)); |
|
4143 |
|
4144 if (f77_exception_encountered) |
|
4145 (*current_liboctave_error_handler) |
|
4146 ("unrecoverable error in zgbtrf"); |
|
4147 else |
|
4148 { |
|
4149 rcond = 0.0; |
|
4150 if (err != 0) |
|
4151 { |
|
4152 err = -2; |
|
4153 |
|
4154 if (sing_handler) |
|
4155 sing_handler (rcond); |
|
4156 else |
|
4157 (*current_liboctave_error_handler) |
|
4158 ("matrix singular to machine precision"); |
|
4159 |
|
4160 } |
|
4161 else |
|
4162 { |
|
4163 char job = 'N'; |
|
4164 int b_nc = b.cols (); |
|
4165 retval = ComplexMatrix (b); |
|
4166 Complex *result = retval.fortran_vec (); |
|
4167 |
|
4168 F77_XFCN (zgbtrs, ZGBTRS, |
|
4169 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4170 nr, n_lower, n_upper, b_nc, tmp_data, |
|
4171 ldm, pipvt, result, b.rows (), err |
|
4172 F77_CHAR_ARG_LEN (1))); |
|
4173 |
|
4174 if (f77_exception_encountered) |
|
4175 { |
|
4176 (*current_liboctave_error_handler) |
|
4177 ("unrecoverable error in dgbtrs"); |
|
4178 } |
|
4179 } |
|
4180 } |
|
4181 } |
|
4182 else if (typ != SparseType::Banded_Hermitian) |
|
4183 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4184 } |
|
4185 |
|
4186 return retval; |
|
4187 } |
|
4188 |
|
4189 SparseComplexMatrix |
|
4190 SparseComplexMatrix::bsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
4191 int& err, double& rcond, |
|
4192 solve_singularity_handler sing_handler) const |
|
4193 { |
|
4194 SparseComplexMatrix retval; |
|
4195 |
|
4196 int nr = rows (); |
|
4197 int nc = cols (); |
|
4198 err = 0; |
|
4199 |
|
4200 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4201 (*current_liboctave_error_handler) |
|
4202 ("matrix dimension mismatch solution of linear equations"); |
|
4203 else |
|
4204 { |
|
4205 // Print spparms("spumoni") info if requested |
|
4206 volatile int typ = mattype.type (); |
|
4207 mattype.info (); |
|
4208 |
|
4209 if (typ == SparseType::Banded_Hermitian) |
|
4210 { |
|
4211 int n_lower = mattype.nlower (); |
|
4212 int ldm = n_lower + 1; |
|
4213 |
|
4214 ComplexMatrix m_band (ldm, nc); |
|
4215 Complex *tmp_data = m_band.fortran_vec (); |
|
4216 |
|
4217 if (! mattype.is_dense ()) |
|
4218 { |
|
4219 int ii = 0; |
|
4220 |
|
4221 for (int j = 0; j < ldm; j++) |
|
4222 for (int i = 0; i < nc; i++) |
|
4223 tmp_data[ii++] = 0.; |
|
4224 } |
|
4225 |
|
4226 for (int j = 0; j < nc; j++) |
|
4227 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4228 { |
|
4229 int ri = ridx (i); |
|
4230 if (ri >= j) |
|
4231 m_band(ri - j, j) = data(i); |
|
4232 } |
|
4233 |
|
4234 char job = 'L'; |
|
4235 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4236 nr, n_lower, tmp_data, ldm, err |
|
4237 F77_CHAR_ARG_LEN (1))); |
|
4238 |
|
4239 if (f77_exception_encountered) |
|
4240 (*current_liboctave_error_handler) |
|
4241 ("unrecoverable error in zpbtrf"); |
|
4242 else |
|
4243 { |
|
4244 rcond = 0.0; |
|
4245 if (err != 0) |
|
4246 { |
|
4247 // Matrix is not positive definite!! Fall through to |
|
4248 // unsymmetric banded solver. |
|
4249 mattype.mark_as_unsymmetric (); |
|
4250 typ = SparseType::Banded; |
|
4251 |
|
4252 err = 0; |
|
4253 } |
|
4254 else |
|
4255 { |
|
4256 rcond = 1.; |
|
4257 int b_nr = b.rows (); |
|
4258 int b_nc = b.cols (); |
|
4259 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); |
|
4260 |
|
4261 // Take a first guess that the number of non-zero terms |
|
4262 // will be as many as in b |
|
4263 volatile int x_nz = b.nnz (); |
|
4264 volatile int ii = 0; |
|
4265 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
4266 |
|
4267 retval.xcidx(0) = 0; |
|
4268 for (volatile int j = 0; j < b_nc; j++) |
|
4269 { |
|
4270 |
|
4271 for (int i = 0; i < b_nr; i++) |
|
4272 Bx[i] = b (i,j); |
|
4273 |
|
4274 F77_XFCN (zpbtrs, ZPBTRS, |
|
4275 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4276 nr, n_lower, 1, tmp_data, |
|
4277 ldm, Bx, b_nr, err |
|
4278 F77_CHAR_ARG_LEN (1))); |
|
4279 |
|
4280 if (f77_exception_encountered) |
|
4281 { |
|
4282 (*current_liboctave_error_handler) |
|
4283 ("unrecoverable error in zpbtrs"); |
|
4284 err = -1; |
|
4285 break; |
|
4286 } |
|
4287 |
|
4288 if (err != 0) |
|
4289 { |
|
4290 (*current_liboctave_error_handler) |
|
4291 ("SparseMatrix::solve solve failed"); |
|
4292 err = -1; |
|
4293 break; |
|
4294 } |
|
4295 |
|
4296 |
|
4297 // Count non-zeros in work vector and adjust |
|
4298 // space in retval if needed |
|
4299 int new_nnz = 0; |
|
4300 for (int i = 0; i < nr; i++) |
|
4301 if (Bx[i] != 0.) |
|
4302 new_nnz++; |
|
4303 |
|
4304 if (ii + new_nnz > x_nz) |
|
4305 { |
|
4306 // Resize the sparse matrix |
|
4307 int sz = new_nnz * (b_nc - j) + x_nz; |
|
4308 retval.change_capacity (sz); |
|
4309 x_nz = sz; |
|
4310 } |
|
4311 |
|
4312 for (int i = 0; i < nr; i++) |
|
4313 if (Bx[i] != 0.) |
|
4314 { |
|
4315 retval.xridx(ii) = i; |
|
4316 retval.xdata(ii++) = Bx[i]; |
|
4317 } |
|
4318 |
|
4319 retval.xcidx(j+1) = ii; |
|
4320 } |
|
4321 |
|
4322 retval.maybe_compress (); |
|
4323 } |
|
4324 } |
|
4325 } |
|
4326 |
|
4327 if (typ == SparseType::Banded) |
|
4328 { |
|
4329 // Create the storage for the banded form of the sparse matrix |
|
4330 int n_upper = mattype.nupper (); |
|
4331 int n_lower = mattype.nlower (); |
|
4332 int ldm = n_upper + 2 * n_lower + 1; |
|
4333 |
|
4334 ComplexMatrix m_band (ldm, nc); |
|
4335 Complex *tmp_data = m_band.fortran_vec (); |
|
4336 |
|
4337 if (! mattype.is_dense ()) |
|
4338 { |
|
4339 int ii = 0; |
|
4340 |
|
4341 for (int j = 0; j < ldm; j++) |
|
4342 for (int i = 0; i < nc; i++) |
|
4343 tmp_data[ii++] = 0.; |
|
4344 } |
|
4345 |
|
4346 for (int j = 0; j < nc; j++) |
|
4347 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4348 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
4349 |
|
4350 Array<int> ipvt (nr); |
|
4351 int *pipvt = ipvt.fortran_vec (); |
|
4352 |
|
4353 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
4354 ldm, pipvt, err)); |
|
4355 |
|
4356 if (f77_exception_encountered) |
|
4357 (*current_liboctave_error_handler) |
|
4358 ("unrecoverable error in xgbtrf"); |
|
4359 else |
|
4360 { |
|
4361 rcond = 0.0; |
|
4362 if (err != 0) |
|
4363 { |
|
4364 err = -2; |
|
4365 |
|
4366 if (sing_handler) |
|
4367 sing_handler (rcond); |
|
4368 else |
|
4369 (*current_liboctave_error_handler) |
|
4370 ("matrix singular to machine precision"); |
|
4371 |
|
4372 } |
|
4373 else |
|
4374 { |
|
4375 char job = 'N'; |
|
4376 volatile int x_nz = b.nnz (); |
|
4377 int b_nc = b.cols (); |
|
4378 retval = SparseComplexMatrix (nr, b_nc, x_nz); |
|
4379 retval.xcidx(0) = 0; |
|
4380 volatile int ii = 0; |
|
4381 |
|
4382 OCTAVE_LOCAL_BUFFER (Complex, Bx, nr); |
|
4383 |
|
4384 for (volatile int j = 0; j < b_nc; j++) |
|
4385 { |
|
4386 for (int i = 0; i < nr; i++) |
|
4387 Bx[i] = 0.; |
|
4388 |
|
4389 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
4390 Bx[b.ridx(i)] = b.data(i); |
|
4391 |
|
4392 F77_XFCN (zgbtrs, ZGBTRS, |
|
4393 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4394 nr, n_lower, n_upper, 1, tmp_data, |
|
4395 ldm, pipvt, Bx, b.rows (), err |
|
4396 F77_CHAR_ARG_LEN (1))); |
|
4397 |
|
4398 if (f77_exception_encountered) |
|
4399 { |
|
4400 (*current_liboctave_error_handler) |
|
4401 ("unrecoverable error in dgbtrs"); |
|
4402 break; |
|
4403 } |
|
4404 |
|
4405 // Count non-zeros in work vector and adjust |
|
4406 // space in retval if needed |
|
4407 int new_nnz = 0; |
|
4408 for (int i = 0; i < nr; i++) |
|
4409 if (Bx[i] != 0.) |
|
4410 new_nnz++; |
|
4411 |
|
4412 if (ii + new_nnz > x_nz) |
|
4413 { |
|
4414 // Resize the sparse matrix |
|
4415 int sz = new_nnz * (b_nc - j) + x_nz; |
|
4416 retval.change_capacity (sz); |
|
4417 x_nz = sz; |
|
4418 } |
|
4419 |
|
4420 for (int i = 0; i < nr; i++) |
|
4421 if (Bx[i] != 0.) |
|
4422 { |
|
4423 retval.xridx(ii) = i; |
|
4424 retval.xdata(ii++) = Bx[i]; |
|
4425 } |
|
4426 retval.xcidx(j+1) = ii; |
|
4427 } |
|
4428 |
|
4429 retval.maybe_compress (); |
|
4430 } |
|
4431 } |
|
4432 } |
|
4433 else if (typ != SparseType::Banded_Hermitian) |
|
4434 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4435 } |
|
4436 |
|
4437 return retval; |
|
4438 } |
|
4439 |
|
4440 void * |
|
4441 SparseComplexMatrix::factorize (int& err, double &rcond, Matrix &Control, |
|
4442 Matrix &Info, |
|
4443 solve_singularity_handler sing_handler) const |
|
4444 { |
|
4445 // The return values |
|
4446 void *Numeric; |
|
4447 err = 0; |
|
4448 |
|
4449 // Setup the control parameters |
|
4450 Control = Matrix (UMFPACK_CONTROL, 1); |
|
4451 double *control = Control.fortran_vec (); |
|
4452 umfpack_zi_defaults (control); |
|
4453 |
|
4454 double tmp = Voctave_sparse_controls.get_key ("spumoni"); |
|
4455 if (!xisnan (tmp)) |
|
4456 Control (UMFPACK_PRL) = tmp; |
|
4457 tmp = Voctave_sparse_controls.get_key ("piv_tol"); |
|
4458 if (!xisnan (tmp)) |
|
4459 { |
|
4460 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
|
4461 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
4462 } |
|
4463 |
|
4464 // Set whether we are allowed to modify Q or not |
|
4465 tmp = Voctave_sparse_controls.get_key ("autoamd"); |
|
4466 if (!xisnan (tmp)) |
|
4467 Control (UMFPACK_FIXQ) = tmp; |
|
4468 |
|
4469 umfpack_zi_report_control (control); |
|
4470 |
|
4471 const int *Ap = cidx (); |
|
4472 const int *Ai = ridx (); |
|
4473 const Complex *Ax = data (); |
|
4474 int nr = rows (); |
|
4475 int nc = cols (); |
|
4476 |
|
4477 umfpack_zi_report_matrix (nr, nc, Ap, Ai, X_CAST (const double *, Ax), |
|
4478 NULL, 1, control); |
|
4479 |
|
4480 void *Symbolic; |
|
4481 Info = Matrix (1, UMFPACK_INFO); |
|
4482 double *info = Info.fortran_vec (); |
|
4483 int status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, |
|
4484 X_CAST (const double *, Ax), |
|
4485 NULL, NULL, &Symbolic, control, info); |
|
4486 |
|
4487 if (status < 0) |
|
4488 { |
|
4489 (*current_liboctave_error_handler) |
|
4490 ("SparseComplexMatrix::solve symbolic factorization failed"); |
|
4491 err = -1; |
|
4492 |
|
4493 umfpack_zi_report_status (control, status); |
|
4494 umfpack_zi_report_info (control, info); |
|
4495 |
|
4496 umfpack_zi_free_symbolic (&Symbolic) ; |
|
4497 } |
|
4498 else |
|
4499 { |
|
4500 umfpack_zi_report_symbolic (Symbolic, control); |
|
4501 |
|
4502 status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), NULL, |
|
4503 Symbolic, &Numeric, control, info) ; |
|
4504 umfpack_zi_free_symbolic (&Symbolic) ; |
|
4505 |
|
4506 #ifdef HAVE_LSSOLVE |
|
4507 rcond = Info (UMFPACK_RCOND); |
|
4508 volatile double rcond_plus_one = rcond + 1.0; |
|
4509 |
|
4510 if (status == UMFPACK_WARNING_singular_matrix || |
|
4511 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
4512 { |
|
4513 umfpack_zi_report_numeric (Numeric, control); |
|
4514 |
|
4515 err = -2; |
|
4516 |
|
4517 if (sing_handler) |
|
4518 sing_handler (rcond); |
|
4519 else |
|
4520 (*current_liboctave_error_handler) |
|
4521 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
4522 rcond); |
|
4523 |
|
4524 } |
|
4525 else |
|
4526 #endif |
|
4527 if (status < 0) |
|
4528 { |
|
4529 (*current_liboctave_error_handler) |
|
4530 ("SparseComplexMatrix::solve numeric factorization failed"); |
|
4531 |
|
4532 umfpack_zi_report_status (control, status); |
|
4533 umfpack_zi_report_info (control, info); |
|
4534 |
|
4535 err = -1; |
|
4536 } |
|
4537 else |
|
4538 { |
|
4539 umfpack_zi_report_numeric (Numeric, control); |
|
4540 } |
|
4541 } |
|
4542 |
|
4543 if (err != 0) |
|
4544 umfpack_zi_free_numeric (&Numeric); |
|
4545 |
|
4546 return Numeric; |
|
4547 } |
|
4548 |
|
4549 ComplexMatrix |
|
4550 SparseComplexMatrix::fsolve (SparseType &mattype, const Matrix& b, int& err, |
|
4551 double& rcond, |
|
4552 solve_singularity_handler sing_handler) const |
|
4553 { |
|
4554 ComplexMatrix retval; |
|
4555 |
|
4556 int nr = rows (); |
|
4557 int nc = cols (); |
|
4558 err = 0; |
|
4559 |
|
4560 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4561 (*current_liboctave_error_handler) |
|
4562 ("matrix dimension mismatch solution of linear equations"); |
|
4563 else |
|
4564 { |
|
4565 // Print spparms("spumoni") info if requested |
|
4566 volatile int typ = mattype.type (); |
|
4567 mattype.info (); |
|
4568 |
|
4569 if (typ == SparseType::Hermitian) |
|
4570 { |
|
4571 // XXX FIXME XXX Write the cholesky solver and only fall |
|
4572 // through if cholesky factorization fails |
|
4573 |
|
4574 (*current_liboctave_warning_handler) |
|
4575 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
4576 |
|
4577 mattype.mark_as_unsymmetric (); |
|
4578 typ = SparseType::Full; |
|
4579 } |
|
4580 |
|
4581 if (typ == SparseType::Full) |
|
4582 { |
|
4583 Matrix Control, Info; |
|
4584 void *Numeric = factorize (err, rcond, Control, Info, |
|
4585 sing_handler); |
|
4586 |
|
4587 if (err == 0) |
|
4588 { |
|
4589 int b_nr = b.rows (); |
|
4590 int b_nc = b.cols (); |
|
4591 int status = 0; |
|
4592 double *control = Control.fortran_vec (); |
|
4593 double *info = Info.fortran_vec (); |
|
4594 const int *Ap = cidx (); |
|
4595 const int *Ai = ridx (); |
|
4596 const Complex *Ax = data (); |
|
4597 const double *Bx = b.fortran_vec (); |
|
4598 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); |
|
4599 for (int i = 0; i < b_nr; i++) |
|
4600 Bz[i] = 0.; |
|
4601 |
|
4602 retval.resize (b_nr, b_nc); |
|
4603 Complex *Xx = retval.fortran_vec (); |
|
4604 |
|
4605 for (int j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) |
|
4606 { |
|
4607 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, |
|
4608 X_CAST (const double *, Ax), |
|
4609 NULL, |
|
4610 X_CAST (double *, &Xx[iidx]), |
|
4611 NULL, |
|
4612 &Bx[iidx], Bz, Numeric, |
|
4613 control, info); |
|
4614 if (status < 0) |
|
4615 { |
|
4616 (*current_liboctave_error_handler) |
|
4617 ("SparseComplexMatrix::solve solve failed"); |
|
4618 |
|
4619 umfpack_zi_report_status (control, status); |
|
4620 |
|
4621 err = -1; |
|
4622 |
|
4623 break; |
|
4624 } |
|
4625 } |
|
4626 |
|
4627 #ifndef HAVE_LSSOLVE |
|
4628 rcond = Info (UMFPACK_RCOND); |
|
4629 volatile double rcond_plus_one = rcond + 1.0; |
|
4630 |
|
4631 if (status == UMFPACK_WARNING_singular_matrix || |
|
4632 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
4633 { |
|
4634 err = -2; |
|
4635 |
|
4636 if (sing_handler) |
|
4637 sing_handler (rcond); |
|
4638 else |
|
4639 (*current_liboctave_error_handler) |
|
4640 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
4641 rcond); |
|
4642 |
|
4643 } |
|
4644 #endif |
|
4645 |
|
4646 umfpack_zi_report_info (control, info); |
|
4647 |
|
4648 umfpack_zi_free_numeric (&Numeric); |
|
4649 } |
|
4650 } |
|
4651 else if (typ != SparseType::Hermitian) |
|
4652 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4653 } |
|
4654 |
|
4655 return retval; |
|
4656 } |
|
4657 |
|
4658 SparseComplexMatrix |
|
4659 SparseComplexMatrix::fsolve (SparseType &mattype, const SparseMatrix& b, |
|
4660 int& err, double& rcond, |
|
4661 solve_singularity_handler sing_handler) const |
|
4662 { |
|
4663 SparseComplexMatrix retval; |
|
4664 |
|
4665 int nr = rows (); |
|
4666 int nc = cols (); |
|
4667 err = 0; |
|
4668 |
|
4669 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4670 (*current_liboctave_error_handler) |
|
4671 ("matrix dimension mismatch solution of linear equations"); |
|
4672 else |
|
4673 { |
|
4674 // Print spparms("spumoni") info if requested |
|
4675 int typ = mattype.type (); |
|
4676 mattype.info (); |
|
4677 |
|
4678 if (typ == SparseType::Hermitian) |
|
4679 { |
|
4680 // XXX FIXME XXX Write the cholesky solver and only fall |
|
4681 // through if cholesky factorization fails |
|
4682 |
|
4683 (*current_liboctave_warning_handler) |
|
4684 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
4685 |
|
4686 mattype.mark_as_unsymmetric (); |
|
4687 typ = SparseType::Full; |
|
4688 } |
|
4689 |
|
4690 if (typ == SparseType::Full) |
|
4691 { |
|
4692 Matrix Control, Info; |
|
4693 void *Numeric = factorize (err, rcond, Control, Info, sing_handler); |
|
4694 |
|
4695 if (err == 0) |
|
4696 { |
|
4697 int b_nr = b.rows (); |
|
4698 int b_nc = b.cols (); |
|
4699 int status = 0; |
|
4700 double *control = Control.fortran_vec (); |
|
4701 double *info = Info.fortran_vec (); |
|
4702 const int *Ap = cidx (); |
|
4703 const int *Ai = ridx (); |
|
4704 const Complex *Ax = data (); |
|
4705 |
|
4706 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
4707 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); |
|
4708 for (int i = 0; i < b_nr; i++) |
|
4709 Bz[i] = 0.; |
|
4710 |
|
4711 // Take a first guess that the number of non-zero terms |
|
4712 // will be as many as in b |
|
4713 int x_nz = b.nnz (); |
|
4714 int ii = 0; |
|
4715 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
4716 |
|
4717 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); |
|
4718 |
|
4719 retval.xcidx(0) = 0; |
|
4720 for (int j = 0; j < b_nc; j++) |
|
4721 { |
|
4722 |
|
4723 for (int i = 0; i < b_nr; i++) |
|
4724 Bx[i] = b.elem (i, j); |
|
4725 |
|
4726 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, |
|
4727 X_CAST (const double *, Ax), |
|
4728 NULL, |
|
4729 X_CAST (double *, Xx), NULL, |
|
4730 Bx, Bz, Numeric, control, |
|
4731 info); |
|
4732 if (status < 0) |
|
4733 { |
|
4734 (*current_liboctave_error_handler) |
|
4735 ("SparseComplexMatrix::solve solve failed"); |
|
4736 |
|
4737 umfpack_zi_report_status (control, status); |
|
4738 |
|
4739 err = -1; |
|
4740 |
|
4741 break; |
|
4742 } |
|
4743 |
|
4744 for (int i = 0; i < b_nr; i++) |
|
4745 { |
|
4746 Complex tmp = Xx[i]; |
|
4747 if (tmp != 0.0) |
|
4748 { |
|
4749 if (ii == x_nz) |
|
4750 { |
|
4751 // Resize the sparse matrix |
|
4752 int sz = x_nz * (b_nc - j) / b_nc; |
|
4753 sz = (sz > 10 ? sz : 10) + x_nz; |
|
4754 retval.change_capacity (sz); |
|
4755 x_nz = sz; |
|
4756 } |
|
4757 retval.xdata(ii) = tmp; |
|
4758 retval.xridx(ii++) = i; |
|
4759 } |
|
4760 } |
|
4761 retval.xcidx(j+1) = ii; |
|
4762 } |
|
4763 |
|
4764 retval.maybe_compress (); |
|
4765 |
|
4766 #ifndef HAVE_LSSOLVE |
|
4767 rcond = Info (UMFPACK_RCOND); |
|
4768 volatile double rcond_plus_one = rcond + 1.0; |
|
4769 |
|
4770 if (status == UMFPACK_WARNING_singular_matrix || |
|
4771 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
4772 { |
|
4773 err = -2; |
|
4774 |
|
4775 if (sing_handler) |
|
4776 sing_handler (rcond); |
|
4777 else |
|
4778 (*current_liboctave_error_handler) |
|
4779 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
4780 rcond); |
|
4781 |
|
4782 } |
|
4783 #endif |
|
4784 |
|
4785 umfpack_zi_report_info (control, info); |
|
4786 |
|
4787 umfpack_zi_free_numeric (&Numeric); |
|
4788 } |
|
4789 } |
|
4790 else if (typ != SparseType::Hermitian) |
|
4791 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4792 } |
|
4793 |
|
4794 return retval; |
|
4795 } |
|
4796 |
|
4797 ComplexMatrix |
|
4798 SparseComplexMatrix::fsolve (SparseType &mattype, const ComplexMatrix& b, |
|
4799 int& err, double& rcond, |
|
4800 solve_singularity_handler sing_handler) const |
|
4801 { |
|
4802 ComplexMatrix retval; |
|
4803 |
|
4804 int nr = rows (); |
|
4805 int nc = cols (); |
|
4806 err = 0; |
|
4807 |
|
4808 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4809 (*current_liboctave_error_handler) |
|
4810 ("matrix dimension mismatch solution of linear equations"); |
|
4811 else |
|
4812 { |
|
4813 // Print spparms("spumoni") info if requested |
|
4814 int typ = mattype.type (); |
|
4815 mattype.info (); |
|
4816 |
|
4817 if (typ == SparseType::Hermitian) |
|
4818 { |
|
4819 // XXX FIXME XXX Write the cholesky solver and only fall |
|
4820 // through if cholesky factorization fails |
|
4821 |
|
4822 (*current_liboctave_warning_handler) |
|
4823 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
4824 |
|
4825 mattype.mark_as_unsymmetric (); |
|
4826 typ = SparseType::Full; |
|
4827 } |
|
4828 |
|
4829 if (typ == SparseType::Full) |
|
4830 { |
|
4831 Matrix Control, Info; |
|
4832 void *Numeric = factorize (err, rcond, Control, Info, sing_handler); |
|
4833 |
|
4834 if (err == 0) |
|
4835 { |
|
4836 int b_nr = b.rows (); |
|
4837 int b_nc = b.cols (); |
|
4838 int status = 0; |
|
4839 double *control = Control.fortran_vec (); |
|
4840 double *info = Info.fortran_vec (); |
|
4841 const int *Ap = cidx (); |
|
4842 const int *Ai = ridx (); |
|
4843 const Complex *Ax = data (); |
|
4844 const Complex *Bx = b.fortran_vec (); |
|
4845 |
|
4846 retval.resize (b_nr, b_nc); |
|
4847 Complex *Xx = retval.fortran_vec (); |
|
4848 |
|
4849 for (int j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) |
|
4850 { |
|
4851 status = |
|
4852 umfpack_zi_solve (UMFPACK_A, Ap, Ai, |
|
4853 X_CAST (const double *, Ax), |
|
4854 NULL, X_CAST (double *, &Xx[iidx]), |
|
4855 NULL, X_CAST (const double *, &Bx[iidx]), |
|
4856 NULL, Numeric, control, info); |
|
4857 |
|
4858 if (status < 0) |
|
4859 { |
|
4860 (*current_liboctave_error_handler) |
|
4861 ("SparseComplexMatrix::solve solve failed"); |
|
4862 |
|
4863 umfpack_zi_report_status (control, status); |
|
4864 |
|
4865 err = -1; |
|
4866 |
|
4867 break; |
|
4868 } |
|
4869 } |
|
4870 |
|
4871 #ifndef HAVE_LSSOLVE |
|
4872 rcond = Info (UMFPACK_RCOND); |
|
4873 volatile double rcond_plus_one = rcond + 1.0; |
|
4874 |
|
4875 if (status == UMFPACK_WARNING_singular_matrix || |
|
4876 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
4877 { |
|
4878 err = -2; |
|
4879 |
|
4880 if (sing_handler) |
|
4881 sing_handler (rcond); |
|
4882 else |
|
4883 (*current_liboctave_error_handler) |
|
4884 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
4885 rcond); |
|
4886 |
|
4887 } |
|
4888 #endif |
|
4889 |
|
4890 umfpack_zi_report_info (control, info); |
|
4891 |
|
4892 umfpack_zi_free_numeric (&Numeric); |
|
4893 } |
|
4894 } |
|
4895 else if (typ != SparseType::Hermitian) |
|
4896 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4897 } |
|
4898 |
|
4899 return retval; |
|
4900 } |
|
4901 |
|
4902 SparseComplexMatrix |
|
4903 SparseComplexMatrix::fsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
4904 int& err, double& rcond, |
|
4905 solve_singularity_handler sing_handler) const |
|
4906 { |
|
4907 SparseComplexMatrix retval; |
|
4908 |
|
4909 int nr = rows (); |
|
4910 int nc = cols (); |
|
4911 err = 0; |
|
4912 |
|
4913 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4914 (*current_liboctave_error_handler) |
|
4915 ("matrix dimension mismatch solution of linear equations"); |
|
4916 else |
|
4917 { |
|
4918 // Print spparms("spumoni") info if requested |
|
4919 int typ = mattype.type (); |
|
4920 mattype.info (); |
|
4921 |
|
4922 if (typ == SparseType::Hermitian) |
|
4923 { |
|
4924 // XXX FIXME XXX Write the cholesky solver and only fall |
|
4925 // through if cholesky factorization fails |
|
4926 |
|
4927 (*current_liboctave_warning_handler) |
|
4928 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
4929 |
|
4930 mattype.mark_as_unsymmetric (); |
|
4931 typ = SparseType::Full; |
|
4932 } |
|
4933 |
|
4934 if (typ == SparseType::Full) |
|
4935 { |
|
4936 Matrix Control, Info; |
|
4937 void *Numeric = factorize (err, rcond, Control, Info, sing_handler); |
|
4938 |
|
4939 if (err == 0) |
|
4940 { |
|
4941 int b_nr = b.rows (); |
|
4942 int b_nc = b.cols (); |
|
4943 int status = 0; |
|
4944 double *control = Control.fortran_vec (); |
|
4945 double *info = Info.fortran_vec (); |
|
4946 const int *Ap = cidx (); |
|
4947 const int *Ai = ridx (); |
|
4948 const Complex *Ax = data (); |
|
4949 |
|
4950 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); |
|
4951 |
|
4952 // Take a first guess that the number of non-zero terms |
|
4953 // will be as many as in b |
|
4954 int x_nz = b.nnz (); |
|
4955 int ii = 0; |
|
4956 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
4957 |
|
4958 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); |
|
4959 |
|
4960 retval.xcidx(0) = 0; |
|
4961 for (int j = 0; j < b_nc; j++) |
|
4962 { |
|
4963 for (int i = 0; i < b_nr; i++) |
|
4964 Bx[i] = b (i,j); |
|
4965 |
|
4966 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, |
|
4967 X_CAST (const double *, Ax), |
|
4968 NULL, X_CAST (double *, Xx), |
|
4969 NULL, X_CAST (double *, Bx), |
|
4970 NULL, Numeric, control, info); |
|
4971 |
|
4972 if (status < 0) |
|
4973 { |
|
4974 (*current_liboctave_error_handler) |
|
4975 ("SparseComplexMatrix::solve solve failed"); |
|
4976 |
|
4977 umfpack_zi_report_status (control, status); |
|
4978 |
|
4979 err = -1; |
|
4980 |
|
4981 break; |
|
4982 } |
|
4983 |
|
4984 for (int i = 0; i < b_nr; i++) |
|
4985 { |
|
4986 Complex tmp = Xx[i]; |
|
4987 if (tmp != 0.0) |
|
4988 { |
|
4989 if (ii == x_nz) |
|
4990 { |
|
4991 // Resize the sparse matrix |
|
4992 int sz = x_nz * (b_nc - j) / b_nc; |
|
4993 sz = (sz > 10 ? sz : 10) + x_nz; |
|
4994 retval.change_capacity (sz); |
|
4995 x_nz = sz; |
|
4996 } |
|
4997 retval.xdata(ii) = tmp; |
|
4998 retval.xridx(ii++) = i; |
|
4999 } |
|
5000 } |
|
5001 retval.xcidx(j+1) = ii; |
|
5002 } |
|
5003 |
|
5004 retval.maybe_compress (); |
|
5005 |
|
5006 #ifndef HAVE_LSSOLVE |
|
5007 rcond = Info (UMFPACK_RCOND); |
|
5008 volatile double rcond_plus_one = rcond + 1.0; |
|
5009 |
|
5010 if (status == UMFPACK_WARNING_singular_matrix || |
|
5011 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
5012 { |
|
5013 err = -2; |
|
5014 |
|
5015 if (sing_handler) |
|
5016 sing_handler (rcond); |
|
5017 else |
|
5018 (*current_liboctave_error_handler) |
|
5019 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", |
|
5020 rcond); |
|
5021 |
|
5022 } |
|
5023 #endif |
|
5024 |
|
5025 umfpack_zi_report_info (control, info); |
|
5026 |
|
5027 umfpack_zi_free_numeric (&Numeric); |
|
5028 } |
|
5029 } |
|
5030 else if (typ != SparseType::Hermitian) |
|
5031 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
5032 } |
|
5033 |
|
5034 return retval; |
|
5035 } |
|
5036 |
|
5037 ComplexMatrix |
|
5038 SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b) const |
|
5039 { |
|
5040 int info; |
|
5041 double rcond; |
|
5042 return solve (mattype, b, info, rcond, 0); |
|
5043 } |
|
5044 |
|
5045 ComplexMatrix |
|
5046 SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, |
|
5047 int& info) const |
|
5048 { |
|
5049 double rcond; |
|
5050 return solve (mattype, b, info, rcond, 0); |
|
5051 } |
|
5052 |
|
5053 ComplexMatrix |
|
5054 SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, int& info, |
|
5055 double& rcond) const |
|
5056 { |
|
5057 return solve (mattype, b, info, rcond, 0); |
|
5058 } |
|
5059 |
|
5060 ComplexMatrix |
|
5061 SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, int& err, |
|
5062 double& rcond, |
|
5063 solve_singularity_handler sing_handler) const |
|
5064 { |
|
5065 int typ = mattype.type (); |
|
5066 |
|
5067 if (typ == SparseType::Unknown) |
|
5068 typ = mattype.type (*this); |
|
5069 |
|
5070 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5071 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5072 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5073 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5074 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5075 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5076 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5077 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5078 else if (typ == SparseType::Tridiagonal || |
|
5079 typ == SparseType::Tridiagonal_Hermitian) |
|
5080 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5081 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5082 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5083 else |
|
5084 { |
|
5085 (*current_liboctave_error_handler) |
|
5086 ("matrix dimension mismatch solution of linear equations"); |
|
5087 return ComplexMatrix (); |
|
5088 } |
|
5089 } |
|
5090 |
|
5091 SparseComplexMatrix |
|
5092 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b) const |
|
5093 { |
|
5094 int info; |
|
5095 double rcond; |
|
5096 return solve (mattype, b, info, rcond, 0); |
|
5097 } |
|
5098 |
|
5099 SparseComplexMatrix |
|
5100 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b, |
|
5101 int& info) const |
|
5102 { |
|
5103 double rcond; |
|
5104 return solve (mattype, b, info, rcond, 0); |
|
5105 } |
|
5106 |
|
5107 SparseComplexMatrix |
|
5108 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b, |
|
5109 int& info, double& rcond) const |
|
5110 { |
|
5111 return solve (mattype, b, info, rcond, 0); |
|
5112 } |
|
5113 |
|
5114 SparseComplexMatrix |
|
5115 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b, |
|
5116 int& err, double& rcond, |
|
5117 solve_singularity_handler sing_handler) const |
|
5118 { |
|
5119 int typ = mattype.type (); |
|
5120 |
|
5121 if (typ == SparseType::Unknown) |
|
5122 typ = mattype.type (*this); |
|
5123 |
|
5124 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5125 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5126 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5127 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5128 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5129 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5130 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5131 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5132 else if (typ == SparseType::Tridiagonal || |
|
5133 typ == SparseType::Tridiagonal_Hermitian) |
|
5134 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5135 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5136 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5137 else |
|
5138 { |
|
5139 (*current_liboctave_error_handler) |
|
5140 ("matrix dimension mismatch solution of linear equations"); |
|
5141 return SparseComplexMatrix (); |
|
5142 } |
|
5143 } |
|
5144 |
|
5145 ComplexMatrix |
|
5146 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b) const |
|
5147 { |
|
5148 int info; |
|
5149 double rcond; |
|
5150 return solve (mattype, b, info, rcond, 0); |
|
5151 } |
|
5152 |
|
5153 ComplexMatrix |
|
5154 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b, |
|
5155 int& info) const |
|
5156 { |
|
5157 double rcond; |
|
5158 return solve (mattype, b, info, rcond, 0); |
|
5159 } |
|
5160 |
|
5161 ComplexMatrix |
|
5162 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b, |
|
5163 int& info, double& rcond) const |
|
5164 { |
|
5165 return solve (mattype, b, info, rcond, 0); |
|
5166 } |
|
5167 |
|
5168 ComplexMatrix |
|
5169 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b, |
|
5170 int& err, double& rcond, |
|
5171 solve_singularity_handler sing_handler) const |
|
5172 { |
|
5173 int typ = mattype.type (); |
|
5174 |
|
5175 if (typ == SparseType::Unknown) |
|
5176 typ = mattype.type (*this); |
|
5177 |
|
5178 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5179 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5180 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5181 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5182 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5183 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5184 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5185 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5186 else if (typ == SparseType::Tridiagonal || |
|
5187 typ == SparseType::Tridiagonal_Hermitian) |
|
5188 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5189 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5190 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5191 else |
|
5192 { |
|
5193 (*current_liboctave_error_handler) |
|
5194 ("matrix dimension mismatch solution of linear equations"); |
|
5195 return ComplexMatrix (); |
|
5196 } |
|
5197 } |
|
5198 |
|
5199 SparseComplexMatrix |
|
5200 SparseComplexMatrix::solve (SparseType &mattype, |
|
5201 const SparseComplexMatrix& b) const |
|
5202 { |
|
5203 int info; |
|
5204 double rcond; |
|
5205 return solve (mattype, b, info, rcond, 0); |
|
5206 } |
|
5207 |
|
5208 SparseComplexMatrix |
|
5209 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, |
|
5210 int& info) const |
|
5211 { |
|
5212 double rcond; |
|
5213 return solve (mattype, b, info, rcond, 0); |
|
5214 } |
|
5215 |
|
5216 SparseComplexMatrix |
|
5217 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, |
|
5218 int& info, double& rcond) const |
|
5219 { |
|
5220 return solve (mattype, b, info, rcond, 0); |
|
5221 } |
|
5222 |
|
5223 SparseComplexMatrix |
|
5224 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, |
|
5225 int& err, double& rcond, |
|
5226 solve_singularity_handler sing_handler) const |
|
5227 { |
|
5228 int typ = mattype.type (); |
|
5229 |
|
5230 if (typ == SparseType::Unknown) |
|
5231 typ = mattype.type (*this); |
|
5232 |
|
5233 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5234 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5235 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5236 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5237 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5238 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5239 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5240 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5241 else if (typ == SparseType::Tridiagonal || |
|
5242 typ == SparseType::Tridiagonal_Hermitian) |
|
5243 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5244 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5245 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5246 else |
|
5247 { |
|
5248 (*current_liboctave_error_handler) |
|
5249 ("matrix dimension mismatch solution of linear equations"); |
|
5250 return SparseComplexMatrix (); |
|
5251 } |
|
5252 } |
|
5253 |
|
5254 ComplexColumnVector |
|
5255 SparseComplexMatrix::solve (SparseType &mattype, const ColumnVector& b) const |
|
5256 { |
|
5257 int info; double rcond; |
|
5258 return solve (mattype, b, info, rcond); |
|
5259 } |
|
5260 |
|
5261 ComplexColumnVector |
|
5262 SparseComplexMatrix::solve (SparseType &mattype, const ColumnVector& b, |
|
5263 int& info) const |
|
5264 { |
|
5265 double rcond; |
|
5266 return solve (mattype, b, info, rcond); |
|
5267 } |
|
5268 |
|
5269 ComplexColumnVector |
|
5270 SparseComplexMatrix::solve (SparseType &mattype, const ColumnVector& b, |
|
5271 int& info, double& rcond) const |
|
5272 { |
|
5273 return solve (mattype, b, info, rcond, 0); |
|
5274 } |
|
5275 |
|
5276 ComplexColumnVector |
|
5277 SparseComplexMatrix::solve (SparseType &mattype, const ColumnVector& b, |
|
5278 int& info, double& rcond, |
|
5279 solve_singularity_handler sing_handler) const |
|
5280 { |
|
5281 Matrix tmp (b); |
|
5282 return solve (mattype, tmp, info, rcond, sing_handler).column (0); |
|
5283 } |
|
5284 |
|
5285 ComplexColumnVector |
|
5286 SparseComplexMatrix::solve (SparseType &mattype, |
|
5287 const ComplexColumnVector& b) const |
|
5288 { |
|
5289 int info; |
|
5290 double rcond; |
|
5291 return solve (mattype, b, info, rcond, 0); |
|
5292 } |
|
5293 |
|
5294 ComplexColumnVector |
|
5295 SparseComplexMatrix::solve (SparseType &mattype, const ComplexColumnVector& b, |
|
5296 int& info) const |
|
5297 { |
|
5298 double rcond; |
|
5299 return solve (mattype, b, info, rcond, 0); |
|
5300 } |
|
5301 |
|
5302 ComplexColumnVector |
|
5303 SparseComplexMatrix::solve (SparseType &mattype, const ComplexColumnVector& b, |
|
5304 int& info, double& rcond) const |
|
5305 { |
|
5306 return solve (mattype, b, info, rcond, 0); |
|
5307 } |
|
5308 |
|
5309 ComplexColumnVector |
|
5310 SparseComplexMatrix::solve (SparseType &mattype, const ComplexColumnVector& b, |
|
5311 int& info, double& rcond, |
|
5312 solve_singularity_handler sing_handler) const |
|
5313 { |
|
5314 ComplexMatrix tmp (b); |
|
5315 return solve (mattype, tmp, info, rcond, sing_handler).column (0); |
|
5316 } |
|
5317 |
|
5318 ComplexMatrix |
|
5319 SparseComplexMatrix::solve (const Matrix& b) const |
|
5320 { |
|
5321 int info; |
|
5322 double rcond; |
|
5323 return solve (b, info, rcond, 0); |
|
5324 } |
|
5325 |
|
5326 ComplexMatrix |
|
5327 SparseComplexMatrix::solve (const Matrix& b, int& info) const |
|
5328 { |
|
5329 double rcond; |
|
5330 return solve (b, info, rcond, 0); |
|
5331 } |
|
5332 |
|
5333 ComplexMatrix |
|
5334 SparseComplexMatrix::solve (const Matrix& b, int& info, |
|
5335 double& rcond) const |
|
5336 { |
|
5337 return solve (b, info, rcond, 0); |
|
5338 } |
|
5339 |
|
5340 ComplexMatrix |
|
5341 SparseComplexMatrix::solve (const Matrix& b, int& err, |
|
5342 double& rcond, |
|
5343 solve_singularity_handler sing_handler) const |
|
5344 { |
|
5345 SparseType mattype (*this); |
|
5346 return solve (mattype, b, err, rcond, sing_handler); |
|
5347 } |
|
5348 |
|
5349 SparseComplexMatrix |
|
5350 SparseComplexMatrix::solve (const SparseMatrix& b) const |
|
5351 { |
|
5352 int info; |
|
5353 double rcond; |
|
5354 return solve (b, info, rcond, 0); |
|
5355 } |
|
5356 |
|
5357 SparseComplexMatrix |
|
5358 SparseComplexMatrix::solve (const SparseMatrix& b, |
|
5359 int& info) const |
|
5360 { |
|
5361 double rcond; |
|
5362 return solve (b, info, rcond, 0); |
|
5363 } |
|
5364 |
|
5365 SparseComplexMatrix |
|
5366 SparseComplexMatrix::solve (const SparseMatrix& b, |
|
5367 int& info, double& rcond) const |
|
5368 { |
|
5369 return solve (b, info, rcond, 0); |
|
5370 } |
|
5371 |
|
5372 SparseComplexMatrix |
|
5373 SparseComplexMatrix::solve (const SparseMatrix& b, |
|
5374 int& err, double& rcond, |
|
5375 solve_singularity_handler sing_handler) const |
|
5376 { |
|
5377 SparseType mattype (*this); |
|
5378 return solve (mattype, b, err, rcond, sing_handler); |
|
5379 } |
|
5380 |
|
5381 ComplexMatrix |
|
5382 SparseComplexMatrix::solve (const ComplexMatrix& b, |
|
5383 int& info) const |
|
5384 { |
|
5385 double rcond; |
|
5386 return solve (b, info, rcond, 0); |
|
5387 } |
|
5388 |
|
5389 ComplexMatrix |
|
5390 SparseComplexMatrix::solve (const ComplexMatrix& b, |
|
5391 int& info, double& rcond) const |
|
5392 { |
|
5393 return solve (b, info, rcond, 0); |
|
5394 } |
|
5395 |
|
5396 ComplexMatrix |
|
5397 SparseComplexMatrix::solve (const ComplexMatrix& b, |
|
5398 int& err, double& rcond, |
|
5399 solve_singularity_handler sing_handler) const |
|
5400 { |
|
5401 SparseType mattype (*this); |
|
5402 return solve (mattype, b, err, rcond, sing_handler); |
|
5403 } |
|
5404 |
|
5405 SparseComplexMatrix |
|
5406 SparseComplexMatrix::solve (const SparseComplexMatrix& b) const |
|
5407 { |
|
5408 int info; |
|
5409 double rcond; |
|
5410 return solve (b, info, rcond, 0); |
|
5411 } |
|
5412 |
|
5413 SparseComplexMatrix |
|
5414 SparseComplexMatrix::solve (const SparseComplexMatrix& b, |
|
5415 int& info) const |
|
5416 { |
|
5417 double rcond; |
|
5418 return solve (b, info, rcond, 0); |
|
5419 } |
|
5420 |
|
5421 SparseComplexMatrix |
|
5422 SparseComplexMatrix::solve (const SparseComplexMatrix& b, |
|
5423 int& info, double& rcond) const |
|
5424 { |
|
5425 return solve (b, info, rcond, 0); |
|
5426 } |
|
5427 |
|
5428 SparseComplexMatrix |
|
5429 SparseComplexMatrix::solve (const SparseComplexMatrix& b, |
|
5430 int& err, double& rcond, |
|
5431 solve_singularity_handler sing_handler) const |
|
5432 { |
|
5433 SparseType mattype (*this); |
|
5434 return solve (mattype, b, err, rcond, sing_handler); |
|
5435 } |
|
5436 |
|
5437 ComplexColumnVector |
|
5438 SparseComplexMatrix::solve (const ColumnVector& b) const |
|
5439 { |
|
5440 int info; double rcond; |
|
5441 return solve (b, info, rcond); |
|
5442 } |
|
5443 |
|
5444 ComplexColumnVector |
|
5445 SparseComplexMatrix::solve (const ColumnVector& b, int& info) const |
|
5446 { |
|
5447 double rcond; |
|
5448 return solve (b, info, rcond); |
|
5449 } |
|
5450 |
|
5451 ComplexColumnVector |
|
5452 SparseComplexMatrix::solve (const ColumnVector& b, int& info, |
|
5453 double& rcond) const |
|
5454 { |
|
5455 return solve (b, info, rcond, 0); |
|
5456 } |
|
5457 |
|
5458 ComplexColumnVector |
|
5459 SparseComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond, |
|
5460 solve_singularity_handler sing_handler) const |
|
5461 { |
|
5462 Matrix tmp (b); |
|
5463 return solve (tmp, info, rcond, sing_handler).column (0); |
|
5464 } |
|
5465 |
|
5466 ComplexColumnVector |
|
5467 SparseComplexMatrix::solve (const ComplexColumnVector& b) const |
|
5468 { |
|
5469 int info; |
|
5470 double rcond; |
|
5471 return solve (b, info, rcond, 0); |
|
5472 } |
|
5473 |
|
5474 ComplexColumnVector |
|
5475 SparseComplexMatrix::solve (const ComplexColumnVector& b, int& info) const |
|
5476 { |
|
5477 double rcond; |
|
5478 return solve (b, info, rcond, 0); |
|
5479 } |
|
5480 |
|
5481 ComplexColumnVector |
|
5482 SparseComplexMatrix::solve (const ComplexColumnVector& b, int& info, |
|
5483 double& rcond) const |
|
5484 { |
|
5485 return solve (b, info, rcond, 0); |
|
5486 } |
|
5487 |
|
5488 ComplexColumnVector |
|
5489 SparseComplexMatrix::solve (const ComplexColumnVector& b, int& info, |
|
5490 double& rcond, |
|
5491 solve_singularity_handler sing_handler) const |
|
5492 { |
|
5493 ComplexMatrix tmp (b); |
|
5494 return solve (tmp, info, rcond, sing_handler).column (0); |
|
5495 } |
|
5496 |
|
5497 ComplexMatrix |
|
5498 SparseComplexMatrix::lssolve (const Matrix& b) const |
|
5499 { |
|
5500 int info; |
|
5501 int rank; |
|
5502 return lssolve (b, info, rank); |
|
5503 } |
|
5504 |
|
5505 ComplexMatrix |
|
5506 SparseComplexMatrix::lssolve (const Matrix& b, int& info) const |
|
5507 { |
|
5508 int rank; |
|
5509 return lssolve (b, info, rank); |
|
5510 } |
|
5511 |
|
5512 ComplexMatrix |
|
5513 SparseComplexMatrix::lssolve (const Matrix& b, int& info, int& rank) const |
|
5514 { |
|
5515 info = -1; |
|
5516 (*current_liboctave_error_handler) |
|
5517 ("SparseComplexMatrix::lssolve not implemented yet"); |
|
5518 return ComplexMatrix (); |
|
5519 } |
|
5520 |
|
5521 SparseComplexMatrix |
|
5522 SparseComplexMatrix::lssolve (const SparseMatrix& b) const |
|
5523 { |
|
5524 int info; |
|
5525 int rank; |
|
5526 return lssolve (b, info, rank); |
|
5527 } |
|
5528 |
|
5529 SparseComplexMatrix |
|
5530 SparseComplexMatrix::lssolve (const SparseMatrix& b, int& info) const |
|
5531 { |
|
5532 int rank; |
|
5533 return lssolve (b, info, rank); |
|
5534 } |
|
5535 |
|
5536 SparseComplexMatrix |
|
5537 SparseComplexMatrix::lssolve (const SparseMatrix& b, int& info, |
|
5538 int& rank) const |
|
5539 { |
|
5540 info = -1; |
|
5541 (*current_liboctave_error_handler) |
|
5542 ("SparseComplexMatrix::lssolve not implemented yet"); |
|
5543 return SparseComplexMatrix (); |
|
5544 } |
|
5545 |
|
5546 ComplexMatrix |
|
5547 SparseComplexMatrix::lssolve (const ComplexMatrix& b) const |
|
5548 { |
|
5549 int info; |
|
5550 int rank; |
|
5551 return lssolve (b, info, rank); |
|
5552 } |
|
5553 |
|
5554 ComplexMatrix |
|
5555 SparseComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const |
|
5556 { |
|
5557 int rank; |
|
5558 return lssolve (b, info, rank); |
|
5559 } |
|
5560 |
|
5561 ComplexMatrix |
|
5562 SparseComplexMatrix::lssolve (const ComplexMatrix& b, int& info, |
|
5563 int& rank) const |
|
5564 { |
|
5565 info = -1; |
|
5566 (*current_liboctave_error_handler) |
|
5567 ("SparseComplexMatrix::lssolve not implemented yet"); |
|
5568 return ComplexMatrix (); |
|
5569 } |
|
5570 |
|
5571 SparseComplexMatrix |
|
5572 SparseComplexMatrix::lssolve (const SparseComplexMatrix& b) const |
|
5573 { |
|
5574 int info; |
|
5575 int rank; |
|
5576 return lssolve (b, info, rank); |
|
5577 } |
|
5578 |
|
5579 SparseComplexMatrix |
|
5580 SparseComplexMatrix::lssolve (const SparseComplexMatrix& b, int& info) const |
|
5581 { |
|
5582 int rank; |
|
5583 return lssolve (b, info, rank); |
|
5584 } |
|
5585 |
|
5586 SparseComplexMatrix |
|
5587 SparseComplexMatrix::lssolve (const SparseComplexMatrix& b, int& info, |
|
5588 int& rank) const |
|
5589 { |
|
5590 info = -1; |
|
5591 (*current_liboctave_error_handler) |
|
5592 ("SparseComplexMatrix::lssolve not implemented yet"); |
|
5593 return SparseComplexMatrix (); |
|
5594 } |
|
5595 |
|
5596 ComplexColumnVector |
|
5597 SparseComplexMatrix::lssolve (const ColumnVector& b) const |
|
5598 { |
|
5599 int info; |
|
5600 int rank; |
|
5601 return lssolve (b, info, rank); |
|
5602 } |
|
5603 |
|
5604 ComplexColumnVector |
|
5605 SparseComplexMatrix::lssolve (const ColumnVector& b, int& info) const |
|
5606 { |
|
5607 int rank; |
|
5608 return lssolve (b, info, rank); |
|
5609 } |
|
5610 |
|
5611 ComplexColumnVector |
|
5612 SparseComplexMatrix::lssolve (const ColumnVector& b, int& info, int& rank) const |
|
5613 { |
|
5614 info = -1; |
|
5615 (*current_liboctave_error_handler) |
|
5616 ("SparseComplexMatrix::lssolve not implemented yet"); |
|
5617 return ComplexColumnVector (); |
|
5618 } |
|
5619 |
|
5620 ComplexColumnVector |
|
5621 SparseComplexMatrix::lssolve (const ComplexColumnVector& b) const |
|
5622 { |
|
5623 int info; |
|
5624 int rank; |
|
5625 return lssolve (b, info, rank); |
|
5626 } |
|
5627 |
|
5628 ComplexColumnVector |
|
5629 SparseComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const |
|
5630 { |
|
5631 int rank; |
|
5632 return lssolve (b, info, rank); |
|
5633 } |
|
5634 |
|
5635 ComplexColumnVector |
|
5636 SparseComplexMatrix::lssolve (const ComplexColumnVector& b, int& info, |
|
5637 int& rank) const |
|
5638 { |
|
5639 info = -1; |
|
5640 (*current_liboctave_error_handler) |
|
5641 ("SparseComplexMatrix::lssolve not implemented yet"); |
|
5642 return ComplexColumnVector (); |
|
5643 } |
|
5644 |
|
5645 // unary operations |
|
5646 SparseBoolMatrix |
|
5647 SparseComplexMatrix::operator ! (void) const |
|
5648 { |
|
5649 int nr = rows (); |
|
5650 int nc = cols (); |
|
5651 int nz1 = nnz (); |
|
5652 int nz2 = nr*nc - nz1; |
|
5653 |
|
5654 SparseBoolMatrix r (nr, nc, nz2); |
|
5655 |
|
5656 int ii = 0; |
|
5657 int jj = 0; |
|
5658 r.cidx (0) = 0; |
|
5659 for (int i = 0; i < nc; i++) |
|
5660 { |
|
5661 for (int j = 0; j < nr; j++) |
|
5662 { |
|
5663 if (jj < cidx(i+1) && ridx(jj) == j) |
|
5664 jj++; |
|
5665 else |
|
5666 { |
|
5667 r.data(ii) = true; |
|
5668 r.ridx(ii++) = j; |
|
5669 } |
|
5670 } |
|
5671 r.cidx (i+1) = ii; |
|
5672 } |
|
5673 |
|
5674 return r; |
|
5675 } |
|
5676 |
|
5677 SparseComplexMatrix |
|
5678 SparseComplexMatrix::squeeze (void) const |
|
5679 { |
|
5680 return MSparse<Complex>::squeeze (); |
|
5681 } |
|
5682 |
|
5683 SparseComplexMatrix |
|
5684 SparseComplexMatrix::index (idx_vector& i, int resize_ok) const |
|
5685 { |
|
5686 return MSparse<Complex>::index (i, resize_ok); |
|
5687 } |
|
5688 |
|
5689 SparseComplexMatrix |
|
5690 SparseComplexMatrix::index (idx_vector& i, idx_vector& j, int resize_ok) const |
|
5691 { |
|
5692 return MSparse<Complex>::index (i, j, resize_ok); |
|
5693 } |
|
5694 |
|
5695 SparseComplexMatrix |
|
5696 SparseComplexMatrix::index (Array<idx_vector>& ra_idx, int resize_ok) const |
|
5697 { |
|
5698 return MSparse<Complex>::index (ra_idx, resize_ok); |
|
5699 } |
|
5700 SparseComplexMatrix |
|
5701 SparseComplexMatrix::reshape (const dim_vector& new_dims) const |
|
5702 { |
|
5703 return MSparse<Complex>::reshape (new_dims); |
|
5704 } |
|
5705 |
|
5706 SparseComplexMatrix |
|
5707 SparseComplexMatrix::permute (const Array<int>& vec, bool inv) const |
|
5708 { |
|
5709 return MSparse<Complex>::permute (vec, inv); |
|
5710 } |
|
5711 |
|
5712 SparseComplexMatrix |
|
5713 SparseComplexMatrix::ipermute (const Array<int>& vec) const |
|
5714 { |
|
5715 return MSparse<Complex>::ipermute (vec); |
|
5716 } |
|
5717 |
|
5718 // other operations |
|
5719 |
|
5720 SparseComplexMatrix |
|
5721 SparseComplexMatrix::map (c_c_Mapper f) const |
|
5722 { |
|
5723 int nr = rows (); |
|
5724 int nc = cols (); |
|
5725 int nz = nnz (); |
|
5726 bool f_zero = (f(0.0) == 0.0); |
|
5727 |
|
5728 // Count number of non-zero elements |
|
5729 int nel = (f_zero ? 0 : nr*nc - nz); |
|
5730 for (int i = 0; i < nz; i++) |
|
5731 if (f (data(i)) != 0.0) |
|
5732 nel++; |
|
5733 |
|
5734 SparseComplexMatrix retval (nr, nc, nel); |
|
5735 |
|
5736 if (f_zero) |
|
5737 { |
|
5738 int ii = 0; |
|
5739 for (int j = 0; j < nc; j++) |
|
5740 { |
|
5741 for (int i = 0; i < nr; i++) |
|
5742 { |
|
5743 Complex tmp = f (elem (i, j)); |
|
5744 if (tmp != 0.0) |
|
5745 { |
|
5746 retval.data(ii) = tmp; |
|
5747 retval.ridx(ii++) = i; |
|
5748 } |
|
5749 } |
|
5750 retval.cidx(j+1) = ii; |
|
5751 } |
|
5752 } |
|
5753 else |
|
5754 { |
|
5755 int ii = 0; |
|
5756 for (int j = 0; j < nc; j++) |
|
5757 { |
|
5758 for (int i = cidx(j); i < cidx(j+1); i++) |
|
5759 { |
|
5760 retval.data(ii) = f (elem(i)); |
|
5761 retval.ridx(ii++) = ridx(i); |
|
5762 } |
|
5763 retval.cidx(j+1) = ii; |
|
5764 } |
|
5765 } |
|
5766 |
|
5767 return retval; |
|
5768 } |
|
5769 |
|
5770 SparseMatrix |
|
5771 SparseComplexMatrix::map (d_c_Mapper f) const |
|
5772 { |
|
5773 int nr = rows (); |
|
5774 int nc = cols (); |
|
5775 int nz = nnz (); |
|
5776 bool f_zero = (f(0.0) == 0.0); |
|
5777 |
|
5778 // Count number of non-zero elements |
|
5779 int nel = (f_zero ? 0 : nr*nc - nz); |
|
5780 for (int i = 0; i < nz; i++) |
|
5781 if (f (data(i)) != 0.0) |
|
5782 nel++; |
|
5783 |
|
5784 SparseMatrix retval (nr, nc, nel); |
|
5785 |
|
5786 if (f_zero) |
|
5787 { |
|
5788 int ii = 0; |
|
5789 for (int j = 0; j < nc; j++) |
|
5790 { |
|
5791 for (int i = 0; i < nr; i++) |
|
5792 { |
|
5793 double tmp = f (elem (i, j)); |
|
5794 if (tmp != 0.0) |
|
5795 { |
|
5796 retval.data(ii) = tmp; |
|
5797 retval.ridx(ii++) = i; |
|
5798 } |
|
5799 } |
|
5800 retval.cidx(j+1) = ii; |
|
5801 } |
|
5802 } |
|
5803 else |
|
5804 { |
|
5805 int ii = 0; |
|
5806 for (int j = 0; j < nc; j++) |
|
5807 { |
|
5808 for (int i = cidx(j); i < cidx(j+1); i++) |
|
5809 { |
|
5810 retval.data(ii) = f (elem(i)); |
|
5811 retval.ridx(ii++) = ridx(i); |
|
5812 } |
|
5813 retval.cidx(j+1) = ii; |
|
5814 } |
|
5815 } |
|
5816 |
|
5817 return retval; |
|
5818 } |
|
5819 |
|
5820 SparseBoolMatrix |
|
5821 SparseComplexMatrix::map (b_c_Mapper f) const |
|
5822 { |
|
5823 int nr = rows (); |
|
5824 int nc = cols (); |
|
5825 int nz = nnz (); |
|
5826 bool f_zero = f(0.0); |
|
5827 |
|
5828 // Count number of non-zero elements |
|
5829 int nel = (f_zero ? 0 : nr*nc - nz); |
|
5830 for (int i = 0; i < nz; i++) |
|
5831 if (f (data(i)) != 0.0) |
|
5832 nel++; |
|
5833 |
|
5834 SparseBoolMatrix retval (nr, nc, nel); |
|
5835 |
|
5836 if (f_zero) |
|
5837 { |
|
5838 int ii = 0; |
|
5839 for (int j = 0; j < nc; j++) |
|
5840 { |
|
5841 for (int i = 0; i < nr; i++) |
|
5842 { |
|
5843 bool tmp = f (elem (i, j)); |
|
5844 if (tmp) |
|
5845 { |
|
5846 retval.data(ii) = tmp; |
|
5847 retval.ridx(ii++) = i; |
|
5848 } |
|
5849 } |
|
5850 retval.cidx(j+1) = ii; |
|
5851 } |
|
5852 } |
|
5853 else |
|
5854 { |
|
5855 int ii = 0; |
|
5856 for (int j = 0; j < nc; j++) |
|
5857 { |
|
5858 for (int i = cidx(j); i < cidx(j+1); i++) |
|
5859 { |
|
5860 retval.data(ii) = f (elem(i)); |
|
5861 retval.ridx(ii++) = ridx(i); |
|
5862 } |
|
5863 retval.cidx(j+1) = ii; |
|
5864 } |
|
5865 } |
|
5866 |
|
5867 return retval; |
|
5868 } |
|
5869 |
|
5870 SparseComplexMatrix& |
|
5871 SparseComplexMatrix::apply (c_c_Mapper f) |
|
5872 { |
|
5873 *this = map (f); |
|
5874 return *this; |
|
5875 } |
|
5876 |
|
5877 bool |
|
5878 SparseComplexMatrix::any_element_is_inf_or_nan (void) const |
|
5879 { |
|
5880 int nel = nnz (); |
|
5881 |
|
5882 for (int i = 0; i < nel; i++) |
|
5883 { |
|
5884 Complex val = data (i); |
|
5885 if (xisinf (val) || xisnan (val)) |
|
5886 return true; |
|
5887 } |
|
5888 |
|
5889 return false; |
|
5890 } |
|
5891 |
|
5892 // Return true if no elements have imaginary components. |
|
5893 |
|
5894 bool |
|
5895 SparseComplexMatrix::all_elements_are_real (void) const |
|
5896 { |
|
5897 int nel = nnz (); |
|
5898 |
|
5899 for (int i = 0; i < nel; i++) |
|
5900 { |
|
5901 double ip = imag (data (i)); |
|
5902 |
|
5903 if (ip != 0.0 || lo_ieee_signbit (ip)) |
|
5904 return false; |
|
5905 } |
|
5906 |
|
5907 return true; |
|
5908 } |
|
5909 |
|
5910 // Return nonzero if any element of CM has a non-integer real or |
|
5911 // imaginary part. Also extract the largest and smallest (real or |
|
5912 // imaginary) values and return them in MAX_VAL and MIN_VAL. |
|
5913 |
|
5914 bool |
|
5915 SparseComplexMatrix::all_integers (double& max_val, double& min_val) const |
|
5916 { |
|
5917 int nel = nnz (); |
|
5918 |
|
5919 if (nel == 0) |
|
5920 return false; |
|
5921 |
|
5922 max_val = real(data (0)); |
|
5923 min_val = real(data (0)); |
|
5924 |
|
5925 for (int i = 0; i < nel; i++) |
|
5926 { |
|
5927 Complex val = data (i); |
|
5928 |
|
5929 double r_val = real (val); |
|
5930 double i_val = imag (val); |
|
5931 |
|
5932 if (r_val > max_val) |
|
5933 max_val = r_val; |
|
5934 |
|
5935 if (i_val > max_val) |
|
5936 max_val = i_val; |
|
5937 |
|
5938 if (r_val < min_val) |
|
5939 min_val = r_val; |
|
5940 |
|
5941 if (i_val < min_val) |
|
5942 min_val = i_val; |
|
5943 |
|
5944 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val) |
|
5945 return false; |
|
5946 } |
|
5947 |
|
5948 return true; |
|
5949 } |
|
5950 |
|
5951 bool |
|
5952 SparseComplexMatrix::too_large_for_float (void) const |
|
5953 { |
|
5954 int nel = nnz (); |
|
5955 |
|
5956 for (int i = 0; i < nel; i++) |
|
5957 { |
|
5958 Complex val = data (i); |
|
5959 |
|
5960 double r_val = real (val); |
|
5961 double i_val = imag (val); |
|
5962 |
|
5963 if (r_val > FLT_MAX |
|
5964 || i_val > FLT_MAX |
|
5965 || r_val < FLT_MIN |
|
5966 || i_val < FLT_MIN) |
|
5967 return true; |
|
5968 } |
|
5969 |
|
5970 return false; |
|
5971 } |
|
5972 |
|
5973 // XXX FIXME XXX Do these really belong here? Maybe they should be |
|
5974 // in a base class? |
|
5975 |
|
5976 SparseBoolMatrix |
|
5977 SparseComplexMatrix::all (int dim) const |
|
5978 { |
|
5979 SPARSE_ALL_OP (dim); |
|
5980 } |
|
5981 |
|
5982 SparseBoolMatrix |
|
5983 SparseComplexMatrix::any (int dim) const |
|
5984 { |
|
5985 SPARSE_ANY_OP (dim); |
|
5986 } |
|
5987 |
|
5988 SparseComplexMatrix |
|
5989 SparseComplexMatrix::cumprod (int dim) const |
|
5990 { |
|
5991 SPARSE_CUMPROD (SparseComplexMatrix, Complex, cumprod); |
|
5992 } |
|
5993 |
|
5994 SparseComplexMatrix |
|
5995 SparseComplexMatrix::cumsum (int dim) const |
|
5996 { |
|
5997 SPARSE_CUMSUM (SparseComplexMatrix, Complex, cumsum); |
|
5998 } |
|
5999 |
|
6000 SparseComplexMatrix |
|
6001 SparseComplexMatrix::prod (int dim) const |
|
6002 { |
|
6003 SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, *=, 1.0, 1.0); |
|
6004 } |
|
6005 |
|
6006 SparseComplexMatrix |
|
6007 SparseComplexMatrix::sum (int dim) const |
|
6008 { |
|
6009 SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, +=, 0.0, 0.0); |
|
6010 } |
|
6011 |
|
6012 SparseComplexMatrix |
|
6013 SparseComplexMatrix::sumsq (int dim) const |
|
6014 { |
|
6015 #define ROW_EXPR \ |
|
6016 Complex d = elem (i, j); \ |
|
6017 tmp [i] += d * conj (d) |
|
6018 |
|
6019 #define COL_EXPR \ |
|
6020 Complex d = elem (i, j); \ |
|
6021 tmp [j] += d * conj (d) |
|
6022 |
|
6023 SPARSE_BASE_REDUCTION_OP (SparseComplexMatrix, Complex, ROW_EXPR, |
|
6024 COL_EXPR, 0.0, 0.0); |
|
6025 |
|
6026 #undef ROW_EXPR |
|
6027 #undef COL_EXPR |
|
6028 } |
|
6029 |
|
6030 SparseMatrix SparseComplexMatrix::abs (void) const |
|
6031 { |
|
6032 int nz = nnz (); |
|
6033 int nc = cols (); |
|
6034 |
|
6035 SparseMatrix retval (rows(), nc, nz); |
|
6036 |
|
6037 for (int i = 0; i < nc + 1; i++) |
|
6038 retval.cidx (i) = cidx (i); |
|
6039 |
|
6040 for (int i = 0; i < nz; i++) |
|
6041 { |
|
6042 retval.data (i) = ::abs (data (i)); |
|
6043 retval.ridx (i) = ridx (i); |
|
6044 } |
|
6045 |
|
6046 return retval; |
|
6047 } |
|
6048 |
|
6049 SparseComplexMatrix |
|
6050 SparseComplexMatrix::diag (int k) const |
|
6051 { |
|
6052 int nnr = rows (); |
|
6053 int nnc = cols (); |
|
6054 |
|
6055 if (k > 0) |
|
6056 nnc -= k; |
|
6057 else if (k < 0) |
|
6058 nnr += k; |
|
6059 |
|
6060 SparseComplexMatrix d; |
|
6061 |
|
6062 if (nnr > 0 && nnc > 0) |
|
6063 { |
|
6064 int ndiag = (nnr < nnc) ? nnr : nnc; |
|
6065 |
|
6066 // Count the number of non-zero elements |
|
6067 int nel = 0; |
|
6068 if (k > 0) |
|
6069 { |
|
6070 for (int i = 0; i < ndiag; i++) |
|
6071 if (elem (i, i+k) != 0.) |
|
6072 nel++; |
|
6073 } |
|
6074 else if ( k < 0) |
|
6075 { |
|
6076 for (int i = 0; i < ndiag; i++) |
|
6077 if (elem (i-k, i) != 0.) |
|
6078 nel++; |
|
6079 } |
|
6080 else |
|
6081 { |
|
6082 for (int i = 0; i < ndiag; i++) |
|
6083 if (elem (i, i) != 0.) |
|
6084 nel++; |
|
6085 } |
|
6086 |
|
6087 d = SparseComplexMatrix (ndiag, 1, nel); |
|
6088 d.xcidx (0) = 0; |
|
6089 d.xcidx (1) = nel; |
|
6090 |
|
6091 int ii = 0; |
|
6092 if (k > 0) |
|
6093 { |
|
6094 for (int i = 0; i < ndiag; i++) |
|
6095 { |
|
6096 Complex tmp = elem (i, i+k); |
|
6097 if (tmp != 0.) |
|
6098 { |
|
6099 d.xdata (ii) = tmp; |
|
6100 d.xridx (ii++) = i; |
|
6101 } |
|
6102 } |
|
6103 } |
|
6104 else if ( k < 0) |
|
6105 { |
|
6106 for (int i = 0; i < ndiag; i++) |
|
6107 { |
|
6108 Complex tmp = elem (i-k, i); |
|
6109 if (tmp != 0.) |
|
6110 { |
|
6111 d.xdata (ii) = tmp; |
|
6112 d.xridx (ii++) = i; |
|
6113 } |
|
6114 } |
|
6115 } |
|
6116 else |
|
6117 { |
|
6118 for (int i = 0; i < ndiag; i++) |
|
6119 { |
|
6120 Complex tmp = elem (i, i); |
|
6121 if (tmp != 0.) |
|
6122 { |
|
6123 d.xdata (ii) = tmp; |
|
6124 d.xridx (ii++) = i; |
|
6125 } |
|
6126 } |
|
6127 } |
|
6128 } |
|
6129 else |
|
6130 (*current_liboctave_error_handler) |
|
6131 ("diag: requested diagonal out of range"); |
|
6132 |
|
6133 return d; |
|
6134 } |
|
6135 |
|
6136 std::ostream& |
|
6137 operator << (std::ostream& os, const SparseComplexMatrix& a) |
|
6138 { |
|
6139 int nc = a.cols (); |
|
6140 |
|
6141 // add one to the printed indices to go from |
|
6142 // zero-based to one-based arrays |
|
6143 for (int j = 0; j < nc; j++) { |
|
6144 OCTAVE_QUIT; |
|
6145 for (int i = a.cidx(j); i < a.cidx(j+1); i++) { |
|
6146 os << a.ridx(i) + 1 << " " << j + 1 << " "; |
|
6147 octave_write_complex (os, a.data(i)); |
|
6148 os << "\n"; |
|
6149 } |
|
6150 } |
|
6151 |
|
6152 return os; |
|
6153 } |
|
6154 |
|
6155 std::istream& |
|
6156 operator >> (std::istream& is, SparseComplexMatrix& a) |
|
6157 { |
|
6158 int nr = a.rows (); |
|
6159 int nc = a.cols (); |
|
6160 int nz = a.nnz (); |
|
6161 |
|
6162 if (nr < 1 || nc < 1) |
|
6163 is.clear (std::ios::badbit); |
|
6164 else |
|
6165 { |
|
6166 int itmp, jtmp, jold = 0; |
|
6167 Complex tmp; |
|
6168 int ii = 0; |
|
6169 |
|
6170 a.cidx (0) = 0; |
|
6171 for (int i = 0; i < nz; i++) |
|
6172 { |
|
6173 is >> itmp; |
|
6174 itmp--; |
|
6175 is >> jtmp; |
|
6176 jtmp--; |
|
6177 tmp = octave_read_complex (is); |
|
6178 |
|
6179 if (is) |
|
6180 { |
|
6181 if (jold != jtmp) |
|
6182 { |
|
6183 for (int j = jold; j < jtmp; j++) |
|
6184 a.cidx(j+1) = ii; |
|
6185 |
|
6186 jold = jtmp; |
|
6187 } |
|
6188 a.data (ii) = tmp; |
|
6189 a.ridx (ii++) = itmp; |
|
6190 } |
|
6191 else |
|
6192 goto done; |
|
6193 } |
|
6194 |
|
6195 for (int j = jold; j < nc; j++) |
|
6196 a.cidx(j+1) = ii; |
|
6197 } |
|
6198 |
|
6199 done: |
|
6200 |
|
6201 return is; |
|
6202 } |
|
6203 |
|
6204 SparseComplexMatrix |
|
6205 operator * (const SparseComplexMatrix& m, const SparseMatrix& a) |
|
6206 { |
|
6207 SparseComplexMatrix tmp (a); |
|
6208 return m * tmp; |
|
6209 } |
|
6210 |
|
6211 SparseComplexMatrix |
|
6212 operator * (const SparseMatrix& m, const SparseComplexMatrix& a) |
|
6213 { |
|
6214 SparseComplexMatrix tmp (m); |
|
6215 return tmp * a; |
|
6216 } |
|
6217 |
|
6218 SparseComplexMatrix |
|
6219 operator * (const SparseComplexMatrix& m, const SparseComplexMatrix& a) |
|
6220 { |
|
6221 #ifdef HAVE_SPARSE_BLAS |
|
6222 // XXX FIXME XXX Isn't there a sparse BLAS ?? |
|
6223 #else |
|
6224 // Use Andy's sparse matrix multiply function |
|
6225 SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex); |
|
6226 #endif |
|
6227 } |
|
6228 |
|
6229 // XXX FIXME XXX -- it would be nice to share code among the min/max |
|
6230 // functions below. |
|
6231 |
|
6232 #define EMPTY_RETURN_CHECK(T) \ |
|
6233 if (nr == 0 || nc == 0) \ |
|
6234 return T (nr, nc); |
|
6235 |
|
6236 SparseComplexMatrix |
|
6237 min (const Complex& c, const SparseComplexMatrix& m) |
|
6238 { |
|
6239 SparseComplexMatrix result; |
|
6240 |
|
6241 int nr = m.rows (); |
|
6242 int nc = m.columns (); |
|
6243 |
|
6244 EMPTY_RETURN_CHECK (SparseComplexMatrix); |
|
6245 |
|
6246 if (abs(c) == 0.) |
|
6247 return SparseComplexMatrix (nr, nc); |
|
6248 else |
|
6249 { |
|
6250 result = SparseComplexMatrix (m); |
|
6251 |
|
6252 for (int j = 0; j < nc; j++) |
|
6253 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6254 result.data(i) = xmin(c, m.data(i)); |
|
6255 } |
|
6256 |
|
6257 return result; |
|
6258 } |
|
6259 |
|
6260 SparseComplexMatrix |
|
6261 min (const SparseComplexMatrix& m, const Complex& c) |
|
6262 { |
|
6263 return min (c, m); |
|
6264 } |
|
6265 |
|
6266 SparseComplexMatrix |
|
6267 min (const SparseComplexMatrix& a, const SparseComplexMatrix& b) |
|
6268 { |
|
6269 SparseComplexMatrix r; |
|
6270 |
|
6271 if ((a.rows() == b.rows()) && (a.cols() == b.cols())) |
|
6272 { |
|
6273 int a_nr = a.rows (); |
|
6274 int a_nc = a.cols (); |
|
6275 |
|
6276 int b_nr = b.rows (); |
|
6277 int b_nc = b.cols (); |
|
6278 |
|
6279 if (a_nr == 0 || b_nc == 0 || a.nnz () == 0 || b.nnz () == 0) |
|
6280 return SparseComplexMatrix (a_nr, a_nc); |
|
6281 |
|
6282 if (a_nr != b_nr || a_nc != b_nc) |
|
6283 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc); |
|
6284 else |
|
6285 { |
|
6286 r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ())); |
|
6287 |
|
6288 int jx = 0; |
|
6289 r.cidx (0) = 0; |
|
6290 for (int i = 0 ; i < a_nc ; i++) |
|
6291 { |
|
6292 int ja = a.cidx(i); |
|
6293 int ja_max = a.cidx(i+1); |
|
6294 bool ja_lt_max= ja < ja_max; |
|
6295 |
|
6296 int jb = b.cidx(i); |
|
6297 int jb_max = b.cidx(i+1); |
|
6298 bool jb_lt_max = jb < jb_max; |
|
6299 |
|
6300 while (ja_lt_max || jb_lt_max ) |
|
6301 { |
|
6302 OCTAVE_QUIT; |
|
6303 if ((! jb_lt_max) || |
|
6304 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) |
|
6305 { |
|
6306 Complex tmp = xmin (a.data(ja), 0.); |
|
6307 if (tmp != 0.) |
|
6308 { |
|
6309 r.ridx(jx) = a.ridx(ja); |
|
6310 r.data(jx) = tmp; |
|
6311 jx++; |
|
6312 } |
|
6313 ja++; |
|
6314 ja_lt_max= ja < ja_max; |
|
6315 } |
|
6316 else if (( !ja_lt_max ) || |
|
6317 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) |
|
6318 { |
|
6319 Complex tmp = xmin (0., b.data(jb)); |
|
6320 if (tmp != 0.) |
|
6321 { |
|
6322 r.ridx(jx) = b.ridx(jb); |
|
6323 r.data(jx) = tmp; |
|
6324 jx++; |
|
6325 } |
|
6326 jb++; |
|
6327 jb_lt_max= jb < jb_max; |
|
6328 } |
|
6329 else |
|
6330 { |
|
6331 Complex tmp = xmin (a.data(ja), b.data(jb)); |
|
6332 if (tmp != 0.) |
|
6333 { |
|
6334 r.data(jx) = tmp; |
|
6335 r.ridx(jx) = a.ridx(ja); |
|
6336 jx++; |
|
6337 } |
|
6338 ja++; |
|
6339 ja_lt_max= ja < ja_max; |
|
6340 jb++; |
|
6341 jb_lt_max= jb < jb_max; |
|
6342 } |
|
6343 } |
|
6344 r.cidx(i+1) = jx; |
|
6345 } |
|
6346 |
|
6347 r.maybe_compress (); |
|
6348 } |
|
6349 } |
|
6350 else |
|
6351 (*current_liboctave_error_handler) ("matrix size mismatch"); |
|
6352 |
|
6353 return r; |
|
6354 } |
|
6355 |
|
6356 SparseComplexMatrix |
|
6357 max (const Complex& c, const SparseComplexMatrix& m) |
|
6358 { |
|
6359 SparseComplexMatrix result; |
|
6360 |
|
6361 int nr = m.rows (); |
|
6362 int nc = m.columns (); |
|
6363 |
|
6364 EMPTY_RETURN_CHECK (SparseComplexMatrix); |
|
6365 |
|
6366 // Count the number of non-zero elements |
|
6367 if (xmax(c, 0.) != 0.) |
|
6368 { |
|
6369 result = SparseComplexMatrix (nr, nc, c); |
|
6370 for (int j = 0; j < nc; j++) |
|
6371 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6372 result.xdata(m.ridx(i) + j * nr) = xmax (c, m.data(i)); |
|
6373 } |
|
6374 else |
|
6375 result = SparseComplexMatrix (m); |
|
6376 |
|
6377 return result; |
|
6378 } |
|
6379 |
|
6380 SparseComplexMatrix |
|
6381 max (const SparseComplexMatrix& m, const Complex& c) |
|
6382 { |
|
6383 return max (c, m); |
|
6384 } |
|
6385 |
|
6386 SparseComplexMatrix |
|
6387 max (const SparseComplexMatrix& a, const SparseComplexMatrix& b) |
|
6388 { |
|
6389 SparseComplexMatrix r; |
|
6390 |
|
6391 if ((a.rows() == b.rows()) && (a.cols() == b.cols())) |
|
6392 { |
|
6393 int a_nr = a.rows (); |
|
6394 int a_nc = a.cols (); |
|
6395 |
|
6396 int b_nr = b.rows (); |
|
6397 int b_nc = b.cols (); |
|
6398 |
|
6399 if (a_nr == 0 || b_nc == 0) |
|
6400 return SparseComplexMatrix (a_nr, a_nc); |
|
6401 if (a.nnz () == 0) |
|
6402 return SparseComplexMatrix (b); |
|
6403 if (b.nnz () == 0) |
|
6404 return SparseComplexMatrix (a); |
|
6405 |
|
6406 if (a_nr != b_nr || a_nc != b_nc) |
|
6407 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc); |
|
6408 else |
|
6409 { |
|
6410 r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ())); |
|
6411 |
|
6412 int jx = 0; |
|
6413 r.cidx (0) = 0; |
|
6414 for (int i = 0 ; i < a_nc ; i++) |
|
6415 { |
|
6416 int ja = a.cidx(i); |
|
6417 int ja_max = a.cidx(i+1); |
|
6418 bool ja_lt_max= ja < ja_max; |
|
6419 |
|
6420 int jb = b.cidx(i); |
|
6421 int jb_max = b.cidx(i+1); |
|
6422 bool jb_lt_max = jb < jb_max; |
|
6423 |
|
6424 while (ja_lt_max || jb_lt_max ) |
|
6425 { |
|
6426 OCTAVE_QUIT; |
|
6427 if ((! jb_lt_max) || |
|
6428 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) |
|
6429 { |
|
6430 Complex tmp = xmax (a.data(ja), 0.); |
|
6431 if (tmp != 0.) |
|
6432 { |
|
6433 r.ridx(jx) = a.ridx(ja); |
|
6434 r.data(jx) = tmp; |
|
6435 jx++; |
|
6436 } |
|
6437 ja++; |
|
6438 ja_lt_max= ja < ja_max; |
|
6439 } |
|
6440 else if (( !ja_lt_max ) || |
|
6441 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) |
|
6442 { |
|
6443 Complex tmp = xmax (0., b.data(jb)); |
|
6444 if (tmp != 0.) |
|
6445 { |
|
6446 r.ridx(jx) = b.ridx(jb); |
|
6447 r.data(jx) = tmp; |
|
6448 jx++; |
|
6449 } |
|
6450 jb++; |
|
6451 jb_lt_max= jb < jb_max; |
|
6452 } |
|
6453 else |
|
6454 { |
|
6455 Complex tmp = xmax (a.data(ja), b.data(jb)); |
|
6456 if (tmp != 0.) |
|
6457 { |
|
6458 r.data(jx) = tmp; |
|
6459 r.ridx(jx) = a.ridx(ja); |
|
6460 jx++; |
|
6461 } |
|
6462 ja++; |
|
6463 ja_lt_max= ja < ja_max; |
|
6464 jb++; |
|
6465 jb_lt_max= jb < jb_max; |
|
6466 } |
|
6467 } |
|
6468 r.cidx(i+1) = jx; |
|
6469 } |
|
6470 |
|
6471 r.maybe_compress (); |
|
6472 } |
|
6473 } |
|
6474 else |
|
6475 (*current_liboctave_error_handler) ("matrix size mismatch"); |
|
6476 |
|
6477 return r; |
|
6478 } |
|
6479 |
|
6480 SPARSE_SMS_CMP_OPS (SparseComplexMatrix, 0.0, real, Complex, |
|
6481 0.0, real) |
|
6482 SPARSE_SMS_BOOL_OPS (SparseComplexMatrix, Complex, 0.0) |
|
6483 |
|
6484 SPARSE_SSM_CMP_OPS (Complex, 0.0, real, SparseComplexMatrix, |
|
6485 0.0, real) |
|
6486 SPARSE_SSM_BOOL_OPS (Complex, SparseComplexMatrix, 0.0) |
|
6487 |
|
6488 SPARSE_SMSM_CMP_OPS (SparseComplexMatrix, 0.0, real, SparseComplexMatrix, |
|
6489 0.0, real) |
|
6490 SPARSE_SMSM_BOOL_OPS (SparseComplexMatrix, SparseComplexMatrix, 0.0) |
|
6491 |
|
6492 /* |
|
6493 ;;; Local Variables: *** |
|
6494 ;;; mode: C++ *** |
|
6495 ;;; End: *** |
|
6496 */ |