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