Mercurial > octave-nkf
comparison liboctave/CSparse.cc @ 5164:57077d0ddc8e
[project @ 2005-02-25 19:55:24 by jwe]
author | jwe |
---|---|
date | Fri, 25 Feb 2005 19:55:28 +0000 |
parents | |
children | dbeafbc0ff64 |
comparison
equal
deleted
inserted
replaced
5163:9f3299378193 | 5164:57077d0ddc8e |
---|---|
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 */ |