Mercurial > hg > octave-nkf
annotate liboctave/CSparse.cc @ 11938:3b5a99b63570 release-3-0-x
backport sorting fixes
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Tue, 24 Feb 2009 07:53:24 +0100 |
| parents | 3fade00a6ac7 |
| children | 402168152bb9 |
| rev | line source |
|---|---|
| 5164 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 2004, 2005, 2006, 2007 David Bateman |
| 7016 | 4 Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 Andy Adler |
| 5 | |
| 6 This file is part of Octave. | |
| 5164 | 7 |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 5164 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 5164 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include <cfloat> | |
| 29 | |
| 30 #include <iostream> | |
| 31 #include <vector> | |
| 32 | |
| 33 #include "quit.h" | |
| 34 #include "lo-ieee.h" | |
| 35 #include "lo-mappers.h" | |
| 36 #include "f77-fcn.h" | |
| 37 #include "dRowVector.h" | |
| 38 | |
| 39 #include "CSparse.h" | |
| 40 #include "boolSparse.h" | |
| 41 #include "dSparse.h" | |
| 42 #include "oct-spparms.h" | |
| 43 #include "SparseCmplxLU.h" | |
| 5451 | 44 #include "oct-sparse.h" |
| 5506 | 45 #include "sparse-util.h" |
| 46 #include "SparseCmplxCHOL.h" | |
| 5610 | 47 #include "SparseCmplxQR.h" |
| 5164 | 48 |
| 5587 | 49 #include "oct-sort.h" |
| 50 | |
| 5681 | 51 // Define whether to use a basic QR solver or one that uses a Dulmange |
| 52 // Mendelsohn factorization to seperate the problem into under-determined, | |
| 53 // well-determined and over-determined parts and solves them seperately | |
| 54 #ifndef USE_QRSOLVE | |
| 55 #include "sparse-dmsolve.cc" | |
| 56 #endif | |
| 57 | |
| 5164 | 58 // Fortran functions we call. |
| 59 extern "C" | |
| 60 { | |
| 61 F77_RET_T | |
| 5275 | 62 F77_FUNC (zgbtrf, ZGBTRF) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, |
| 63 const octave_idx_type&, Complex*, const octave_idx_type&, octave_idx_type*, octave_idx_type&); | |
| 5164 | 64 |
| 65 F77_RET_T | |
| 5275 | 66 F77_FUNC (zgbtrs, ZGBTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 67 const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, | |
| 68 const Complex*, const octave_idx_type&, | |
| 69 const octave_idx_type*, Complex*, const octave_idx_type&, octave_idx_type& | |
| 5164 | 70 F77_CHAR_ARG_LEN_DECL); |
| 71 | |
| 72 F77_RET_T | |
| 5275 | 73 F77_FUNC (zgbcon, ZGBCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 74 const octave_idx_type&, const octave_idx_type&, Complex*, | |
| 75 const octave_idx_type&, const octave_idx_type*, const double&, | |
| 76 double&, Complex*, double*, octave_idx_type& | |
| 5164 | 77 F77_CHAR_ARG_LEN_DECL); |
| 78 | |
| 79 F77_RET_T | |
| 5275 | 80 F77_FUNC (zpbtrf, ZPBTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 81 const octave_idx_type&, Complex*, const octave_idx_type&, octave_idx_type& | |
| 5164 | 82 F77_CHAR_ARG_LEN_DECL); |
| 83 | |
| 84 F77_RET_T | |
| 5275 | 85 F77_FUNC (zpbtrs, ZPBTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 86 const octave_idx_type&, const octave_idx_type&, Complex*, const octave_idx_type&, | |
| 87 Complex*, const octave_idx_type&, octave_idx_type& | |
| 5164 | 88 F77_CHAR_ARG_LEN_DECL); |
| 89 | |
| 90 F77_RET_T | |
| 5275 | 91 F77_FUNC (zpbcon, ZPBCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 92 const octave_idx_type&, Complex*, const octave_idx_type&, | |
| 5681 | 93 const double&, double&, Complex*, double*, octave_idx_type& |
| 5164 | 94 F77_CHAR_ARG_LEN_DECL); |
| 95 | |
| 96 F77_RET_T | |
| 5275 | 97 F77_FUNC (zgttrf, ZGTTRF) (const octave_idx_type&, Complex*, Complex*, Complex*, |
| 98 Complex*, octave_idx_type*, octave_idx_type&); | |
| 5164 | 99 |
| 100 F77_RET_T | |
| 5275 | 101 F77_FUNC (zgttrs, ZGTTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 102 const octave_idx_type&, const Complex*, const Complex*, | |
| 103 const Complex*, const Complex*, const octave_idx_type*, | |
| 104 Complex *, const octave_idx_type&, octave_idx_type& | |
| 5164 | 105 F77_CHAR_ARG_LEN_DECL); |
| 106 | |
| 107 F77_RET_T | |
| 5322 | 108 F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&, double*, Complex*, |
| 5275 | 109 Complex*, const octave_idx_type&, octave_idx_type&); |
| 5164 | 110 |
| 111 F77_RET_T | |
| 5275 | 112 F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&, Complex*, Complex*, |
| 113 Complex*, Complex*, const octave_idx_type&, octave_idx_type&); | |
| 5164 | 114 } |
| 115 | |
| 116 SparseComplexMatrix::SparseComplexMatrix (const SparseMatrix& a) | |
| 5681 | 117 : MSparse<Complex> (a.rows (), a.cols (), a.nnz ()) |
| 5164 | 118 { |
| 5275 | 119 octave_idx_type nc = cols (); |
| 5681 | 120 octave_idx_type nz = a.nnz (); |
| 5275 | 121 |
| 122 for (octave_idx_type i = 0; i < nc + 1; i++) | |
| 5164 | 123 cidx (i) = a.cidx (i); |
| 124 | |
| 5275 | 125 for (octave_idx_type i = 0; i < nz; i++) |
| 5164 | 126 { |
| 5681 | 127 data (i) = Complex (a.data (i)); |
| 5164 | 128 ridx (i) = a.ridx (i); |
| 129 } | |
| 130 } | |
| 131 | |
| 132 SparseComplexMatrix::SparseComplexMatrix (const SparseBoolMatrix& a) | |
| 5681 | 133 : MSparse<Complex> (a.rows (), a.cols (), a.nnz ()) |
| 5164 | 134 { |
| 5275 | 135 octave_idx_type nc = cols (); |
| 5681 | 136 octave_idx_type nz = a.nnz (); |
| 5275 | 137 |
| 138 for (octave_idx_type i = 0; i < nc + 1; i++) | |
| 5164 | 139 cidx (i) = a.cidx (i); |
| 140 | |
| 5275 | 141 for (octave_idx_type i = 0; i < nz; i++) |
| 5164 | 142 { |
| 5681 | 143 data (i) = Complex (a.data (i)); |
| 5164 | 144 ridx (i) = a.ridx (i); |
| 145 } | |
| 146 } | |
| 147 | |
| 148 bool | |
| 149 SparseComplexMatrix::operator == (const SparseComplexMatrix& a) const | |
| 150 { | |
| 5275 | 151 octave_idx_type nr = rows (); |
| 152 octave_idx_type nc = cols (); | |
| 5681 | 153 octave_idx_type nz = nnz (); |
| 5275 | 154 octave_idx_type nr_a = a.rows (); |
| 155 octave_idx_type nc_a = a.cols (); | |
| 5681 | 156 octave_idx_type nz_a = a.nnz (); |
| 5164 | 157 |
| 158 if (nr != nr_a || nc != nc_a || nz != nz_a) | |
| 159 return false; | |
| 160 | |
| 5275 | 161 for (octave_idx_type i = 0; i < nc + 1; i++) |
| 5164 | 162 if (cidx(i) != a.cidx(i)) |
| 163 return false; | |
| 164 | |
| 5275 | 165 for (octave_idx_type i = 0; i < nz; i++) |
| 5164 | 166 if (data(i) != a.data(i) || ridx(i) != a.ridx(i)) |
| 167 return false; | |
| 168 | |
| 169 return true; | |
| 170 } | |
| 171 | |
| 172 bool | |
| 173 SparseComplexMatrix::operator != (const SparseComplexMatrix& a) const | |
| 174 { | |
| 175 return !(*this == a); | |
| 176 } | |
| 177 | |
| 178 bool | |
| 179 SparseComplexMatrix::is_hermitian (void) const | |
| 180 { | |
| 5275 | 181 octave_idx_type nr = rows (); |
| 182 octave_idx_type nc = cols (); | |
| 5164 | 183 |
| 6207 | 184 if (nr == nc && nr > 0) |
| 5164 | 185 { |
| 6207 | 186 for (octave_idx_type j = 0; j < nc; j++) |
| 187 { | |
| 188 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 189 { | |
| 190 octave_idx_type ri = ridx(i); | |
| 191 | |
| 192 if (ri != j) | |
| 193 { | |
| 194 bool found = false; | |
| 195 | |
| 196 for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++) | |
| 197 { | |
| 198 if (ridx(k) == j) | |
| 199 { | |
| 200 if (data(i) == conj(data(k))) | |
| 201 found = true; | |
| 202 break; | |
| 203 } | |
| 204 } | |
| 205 | |
| 206 if (! found) | |
| 207 return false; | |
| 208 } | |
| 209 } | |
| 210 } | |
| 5164 | 211 |
| 212 return true; | |
| 213 } | |
| 214 | |
| 215 return false; | |
| 216 } | |
| 217 | |
| 218 static const Complex Complex_NaN_result (octave_NaN, octave_NaN); | |
| 219 | |
| 220 SparseComplexMatrix | |
| 221 SparseComplexMatrix::max (int dim) const | |
| 222 { | |
| 5275 | 223 Array2<octave_idx_type> dummy_idx; |
| 5164 | 224 return max (dummy_idx, dim); |
| 225 } | |
| 226 | |
| 227 SparseComplexMatrix | |
| 5275 | 228 SparseComplexMatrix::max (Array2<octave_idx_type>& idx_arg, int dim) const |
| 5164 | 229 { |
| 230 SparseComplexMatrix result; | |
| 231 dim_vector dv = dims (); | |
| 232 | |
| 233 if (dv.numel () == 0 || dim > dv.length () || dim < 0) | |
| 234 return result; | |
| 235 | |
| 5275 | 236 octave_idx_type nr = dv(0); |
| 237 octave_idx_type nc = dv(1); | |
| 5164 | 238 |
| 239 if (dim == 0) | |
| 240 { | |
| 241 idx_arg.resize (1, nc); | |
| 5275 | 242 octave_idx_type nel = 0; |
| 243 for (octave_idx_type j = 0; j < nc; j++) | |
| 5164 | 244 { |
| 245 Complex tmp_max; | |
| 246 double abs_max = octave_NaN; | |
| 5275 | 247 octave_idx_type idx_j = 0; |
| 248 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 5164 | 249 { |
| 250 if (ridx(i) != idx_j) | |
| 251 break; | |
| 252 else | |
| 253 idx_j++; | |
| 254 } | |
| 255 | |
| 256 if (idx_j != nr) | |
| 257 { | |
| 258 tmp_max = 0.; | |
| 259 abs_max = 0.; | |
| 260 } | |
| 261 | |
| 5275 | 262 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) |
| 5164 | 263 { |
| 264 Complex tmp = data (i); | |
| 265 | |
| 5389 | 266 if (xisnan (tmp)) |
| 5164 | 267 continue; |
| 268 | |
| 5261 | 269 double abs_tmp = std::abs (tmp); |
| 5164 | 270 |
| 5389 | 271 if (xisnan (abs_max) || abs_tmp > abs_max) |
| 5164 | 272 { |
| 273 idx_j = ridx (i); | |
| 274 tmp_max = tmp; | |
| 275 abs_max = abs_tmp; | |
| 276 } | |
| 277 } | |
| 278 | |
| 5389 | 279 idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_j; |
| 5164 | 280 if (abs_max != 0.) |
| 281 nel++; | |
| 282 } | |
| 283 | |
| 284 result = SparseComplexMatrix (1, nc, nel); | |
| 285 | |
| 5275 | 286 octave_idx_type ii = 0; |
| 5164 | 287 result.xcidx (0) = 0; |
| 5275 | 288 for (octave_idx_type j = 0; j < nc; j++) |
| 5164 | 289 { |
| 290 Complex tmp = elem (idx_arg(j), j); | |
| 291 if (tmp != 0.) | |
| 292 { | |
| 293 result.xdata (ii) = tmp; | |
| 294 result.xridx (ii++) = 0; | |
| 295 } | |
| 296 result.xcidx (j+1) = ii; | |
| 297 } | |
| 298 } | |
| 299 else | |
| 300 { | |
| 301 idx_arg.resize (nr, 1, 0); | |
| 302 | |
| 5275 | 303 for (octave_idx_type i = cidx(0); i < cidx(1); i++) |
| 5164 | 304 idx_arg.elem(ridx(i)) = -1; |
| 305 | |
| 5275 | 306 for (octave_idx_type j = 0; j < nc; j++) |
| 307 for (octave_idx_type i = 0; i < nr; i++) | |
| 5164 | 308 { |
| 309 if (idx_arg.elem(i) != -1) | |
| 310 continue; | |
| 311 bool found = false; | |
| 5275 | 312 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++) |
| 5164 | 313 if (ridx(k) == i) |
| 314 { | |
| 315 found = true; | |
| 316 break; | |
| 317 } | |
| 318 | |
| 319 if (!found) | |
| 320 idx_arg.elem(i) = j; | |
| 321 | |
| 322 } | |
| 323 | |
| 5275 | 324 for (octave_idx_type j = 0; j < nc; j++) |
| 5164 | 325 { |
| 5275 | 326 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) |
| 5164 | 327 { |
| 5275 | 328 octave_idx_type ir = ridx (i); |
| 329 octave_idx_type ix = idx_arg.elem (ir); | |
| 5164 | 330 Complex tmp = data (i); |
| 331 | |
| 5389 | 332 if (xisnan (tmp)) |
| 5164 | 333 continue; |
| 5261 | 334 else if (ix == -1 || std::abs(tmp) > std::abs(elem (ir, ix))) |
| 5164 | 335 idx_arg.elem (ir) = j; |
| 336 } | |
| 337 } | |
| 338 | |
| 5275 | 339 octave_idx_type nel = 0; |
| 340 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 341 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) |
| 342 nel++; | |
| 343 | |
| 344 result = SparseComplexMatrix (nr, 1, nel); | |
| 345 | |
| 5275 | 346 octave_idx_type ii = 0; |
| 5164 | 347 result.xcidx (0) = 0; |
| 348 result.xcidx (1) = nel; | |
| 5275 | 349 for (octave_idx_type j = 0; j < nr; j++) |
| 5164 | 350 { |
| 351 if (idx_arg(j) == -1) | |
| 352 { | |
| 353 idx_arg(j) = 0; | |
| 354 result.xdata (ii) = Complex_NaN_result; | |
| 355 result.xridx (ii++) = j; | |
| 356 } | |
| 357 else | |
| 358 { | |
| 359 Complex tmp = elem (j, idx_arg(j)); | |
| 360 if (tmp != 0.) | |
| 361 { | |
| 362 result.xdata (ii) = tmp; | |
| 363 result.xridx (ii++) = j; | |
| 364 } | |
| 365 } | |
| 366 } | |
| 367 } | |
| 368 | |
| 369 return result; | |
| 370 } | |
| 371 | |
| 372 SparseComplexMatrix | |
| 373 SparseComplexMatrix::min (int dim) const | |
| 374 { | |
| 5275 | 375 Array2<octave_idx_type> dummy_idx; |
| 5164 | 376 return min (dummy_idx, dim); |
| 377 } | |
| 378 | |
| 379 SparseComplexMatrix | |
| 5275 | 380 SparseComplexMatrix::min (Array2<octave_idx_type>& idx_arg, int dim) const |
| 5164 | 381 { |
| 382 SparseComplexMatrix result; | |
| 383 dim_vector dv = dims (); | |
| 384 | |
| 385 if (dv.numel () == 0 || dim > dv.length () || dim < 0) | |
| 386 return result; | |
| 387 | |
| 5275 | 388 octave_idx_type nr = dv(0); |
| 389 octave_idx_type nc = dv(1); | |
| 5164 | 390 |
| 391 if (dim == 0) | |
| 392 { | |
| 393 idx_arg.resize (1, nc); | |
| 5275 | 394 octave_idx_type nel = 0; |
| 395 for (octave_idx_type j = 0; j < nc; j++) | |
| 5164 | 396 { |
| 397 Complex tmp_min; | |
| 398 double abs_min = octave_NaN; | |
| 5275 | 399 octave_idx_type idx_j = 0; |
| 400 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 5164 | 401 { |
| 402 if (ridx(i) != idx_j) | |
| 403 break; | |
| 404 else | |
| 405 idx_j++; | |
| 406 } | |
| 407 | |
| 408 if (idx_j != nr) | |
| 409 { | |
| 410 tmp_min = 0.; | |
| 411 abs_min = 0.; | |
| 412 } | |
| 413 | |
| 5275 | 414 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) |
| 5164 | 415 { |
| 416 Complex tmp = data (i); | |
| 417 | |
| 5389 | 418 if (xisnan (tmp)) |
| 5164 | 419 continue; |
| 420 | |
| 5261 | 421 double abs_tmp = std::abs (tmp); |
| 5164 | 422 |
| 5389 | 423 if (xisnan (abs_min) || abs_tmp < abs_min) |
| 5164 | 424 { |
| 425 idx_j = ridx (i); | |
| 426 tmp_min = tmp; | |
| 427 abs_min = abs_tmp; | |
| 428 } | |
| 429 } | |
| 430 | |
| 5389 | 431 idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_j; |
| 5164 | 432 if (abs_min != 0.) |
| 433 nel++; | |
| 434 } | |
| 435 | |
| 436 result = SparseComplexMatrix (1, nc, nel); | |
| 437 | |
| 5275 | 438 octave_idx_type ii = 0; |
| 5164 | 439 result.xcidx (0) = 0; |
| 5275 | 440 for (octave_idx_type j = 0; j < nc; j++) |
| 5164 | 441 { |
| 442 Complex tmp = elem (idx_arg(j), j); | |
| 443 if (tmp != 0.) | |
| 444 { | |
| 445 result.xdata (ii) = tmp; | |
| 446 result.xridx (ii++) = 0; | |
| 447 } | |
| 448 result.xcidx (j+1) = ii; | |
| 449 } | |
| 450 } | |
| 451 else | |
| 452 { | |
| 453 idx_arg.resize (nr, 1, 0); | |
| 454 | |
| 5275 | 455 for (octave_idx_type i = cidx(0); i < cidx(1); i++) |
| 5164 | 456 idx_arg.elem(ridx(i)) = -1; |
| 457 | |
| 5275 | 458 for (octave_idx_type j = 0; j < nc; j++) |
| 459 for (octave_idx_type i = 0; i < nr; i++) | |
| 5164 | 460 { |
| 461 if (idx_arg.elem(i) != -1) | |
| 462 continue; | |
| 463 bool found = false; | |
| 5275 | 464 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++) |
| 5164 | 465 if (ridx(k) == i) |
| 466 { | |
| 467 found = true; | |
| 468 break; | |
| 469 } | |
| 470 | |
| 471 if (!found) | |
| 472 idx_arg.elem(i) = j; | |
| 473 | |
| 474 } | |
| 475 | |
| 5275 | 476 for (octave_idx_type j = 0; j < nc; j++) |
| 5164 | 477 { |
| 5275 | 478 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) |
| 5164 | 479 { |
| 5275 | 480 octave_idx_type ir = ridx (i); |
| 481 octave_idx_type ix = idx_arg.elem (ir); | |
| 5164 | 482 Complex tmp = data (i); |
| 483 | |
| 5389 | 484 if (xisnan (tmp)) |
| 5164 | 485 continue; |
| 5261 | 486 else if (ix == -1 || std::abs(tmp) < std::abs(elem (ir, ix))) |
| 5164 | 487 idx_arg.elem (ir) = j; |
| 488 } | |
| 489 } | |
| 490 | |
| 5275 | 491 octave_idx_type nel = 0; |
| 492 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 493 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) |
| 494 nel++; | |
| 495 | |
| 496 result = SparseComplexMatrix (nr, 1, nel); | |
| 497 | |
| 5275 | 498 octave_idx_type ii = 0; |
| 5164 | 499 result.xcidx (0) = 0; |
| 500 result.xcidx (1) = nel; | |
| 5275 | 501 for (octave_idx_type j = 0; j < nr; j++) |
| 5164 | 502 { |
| 503 if (idx_arg(j) == -1) | |
| 504 { | |
| 505 idx_arg(j) = 0; | |
| 506 result.xdata (ii) = Complex_NaN_result; | |
| 507 result.xridx (ii++) = j; | |
| 508 } | |
| 509 else | |
| 510 { | |
| 511 Complex tmp = elem (j, idx_arg(j)); | |
| 512 if (tmp != 0.) | |
| 513 { | |
| 514 result.xdata (ii) = tmp; | |
| 515 result.xridx (ii++) = j; | |
| 516 } | |
| 517 } | |
| 518 } | |
| 519 } | |
| 520 | |
| 521 return result; | |
| 522 } | |
| 523 | |
| 524 // destructive insert/delete/reorder operations | |
| 525 | |
| 526 SparseComplexMatrix& | |
| 5275 | 527 SparseComplexMatrix::insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c) |
| 5164 | 528 { |
| 529 SparseComplexMatrix tmp (a); | |
| 6060 | 530 return insert (tmp /*a*/, r, c); |
| 5164 | 531 } |
| 532 | |
| 533 SparseComplexMatrix& | |
| 5275 | 534 SparseComplexMatrix::insert (const SparseComplexMatrix& a, octave_idx_type r, octave_idx_type c) |
| 5164 | 535 { |
| 536 MSparse<Complex>::insert (a, r, c); | |
| 537 return *this; | |
| 538 } | |
| 539 | |
| 6823 | 540 SparseComplexMatrix& |
| 541 SparseComplexMatrix::insert (const SparseMatrix& a, const Array<octave_idx_type>& indx) | |
| 542 { | |
| 543 SparseComplexMatrix tmp (a); | |
| 544 return insert (tmp /*a*/, indx); | |
| 545 } | |
| 546 | |
| 547 SparseComplexMatrix& | |
| 548 SparseComplexMatrix::insert (const SparseComplexMatrix& a, const Array<octave_idx_type>& indx) | |
| 549 { | |
| 550 MSparse<Complex>::insert (a, indx); | |
| 551 return *this; | |
| 552 } | |
| 553 | |
| 5164 | 554 SparseComplexMatrix |
| 555 SparseComplexMatrix::concat (const SparseComplexMatrix& rb, | |
| 5275 | 556 const Array<octave_idx_type>& ra_idx) |
| 5164 | 557 { |
| 558 // Don't use numel to avoid all possiblity of an overflow | |
| 559 if (rb.rows () > 0 && rb.cols () > 0) | |
| 560 insert (rb, ra_idx(0), ra_idx(1)); | |
| 561 return *this; | |
| 562 } | |
| 563 | |
| 564 SparseComplexMatrix | |
| 5275 | 565 SparseComplexMatrix::concat (const SparseMatrix& rb, const Array<octave_idx_type>& ra_idx) |
| 5164 | 566 { |
| 567 SparseComplexMatrix tmp (rb); | |
| 568 if (rb.rows () > 0 && rb.cols () > 0) | |
| 569 insert (tmp, ra_idx(0), ra_idx(1)); | |
| 570 return *this; | |
| 571 } | |
| 572 | |
| 573 ComplexMatrix | |
| 574 SparseComplexMatrix::matrix_value (void) const | |
| 575 { | |
| 5275 | 576 octave_idx_type nr = rows (); |
| 577 octave_idx_type nc = cols (); | |
| 5164 | 578 ComplexMatrix retval (nr, nc, Complex (0.0, 0.0)); |
| 579 | |
| 5275 | 580 for (octave_idx_type j = 0; j < nc; j++) |
| 581 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 5164 | 582 retval.elem (ridx(i), j) = data (i); |
| 583 | |
| 584 return retval; | |
| 585 } | |
| 586 | |
| 587 SparseComplexMatrix | |
| 588 SparseComplexMatrix::hermitian (void) const | |
| 589 { | |
| 5275 | 590 octave_idx_type nr = rows (); |
| 591 octave_idx_type nc = cols (); | |
| 5681 | 592 octave_idx_type nz = nnz (); |
| 5164 | 593 SparseComplexMatrix retval (nc, nr, nz); |
| 594 | |
| 5648 | 595 OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr + 1); |
| 596 for (octave_idx_type i = 0; i < nr; i++) | |
| 597 w[i] = 0; | |
| 598 for (octave_idx_type i = 0; i < nz; i++) | |
| 599 w[ridx(i)]++; | |
| 600 nz = 0; | |
| 601 for (octave_idx_type i = 0; i < nr; i++) | |
| 5164 | 602 { |
| 5648 | 603 retval.xcidx(i) = nz; |
| 604 nz += w[i]; | |
| 605 w[i] = retval.xcidx(i); | |
| 5164 | 606 } |
| 5648 | 607 retval.xcidx(nr) = nz; |
| 608 w[nr] = nz; | |
| 609 | |
| 610 for (octave_idx_type j = 0; j < nc; j++) | |
| 611 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++) | |
| 612 { | |
| 613 octave_idx_type q = w [ridx(k)]++; | |
| 614 retval.xridx (q) = j; | |
| 615 retval.xdata (q) = conj (data (k)); | |
| 616 } | |
| 5164 | 617 |
| 618 return retval; | |
| 619 } | |
| 620 | |
| 621 SparseComplexMatrix | |
| 622 conj (const SparseComplexMatrix& a) | |
| 623 { | |
| 5275 | 624 octave_idx_type nr = a.rows (); |
| 625 octave_idx_type nc = a.cols (); | |
| 5681 | 626 octave_idx_type nz = a.nnz (); |
| 5164 | 627 SparseComplexMatrix retval (nc, nr, nz); |
| 628 | |
| 5275 | 629 for (octave_idx_type i = 0; i < nc + 1; i++) |
| 5164 | 630 retval.cidx (i) = a.cidx (i); |
| 631 | |
| 5275 | 632 for (octave_idx_type i = 0; i < nz; i++) |
| 5164 | 633 { |
| 634 retval.data (i) = conj (a.data (i)); | |
| 635 retval.ridx (i) = a.ridx (i); | |
| 636 } | |
| 637 | |
| 638 return retval; | |
| 639 } | |
| 640 | |
| 641 SparseComplexMatrix | |
| 642 SparseComplexMatrix::inverse (void) const | |
| 643 { | |
| 5275 | 644 octave_idx_type info; |
| 5164 | 645 double rcond; |
| 5785 | 646 MatrixType mattype (*this); |
| 5506 | 647 return inverse (mattype, info, rcond, 0, 0); |
| 648 } | |
| 649 | |
| 650 SparseComplexMatrix | |
| 5785 | 651 SparseComplexMatrix::inverse (MatrixType& mattype) const |
| 5506 | 652 { |
| 653 octave_idx_type info; | |
| 654 double rcond; | |
| 655 return inverse (mattype, info, rcond, 0, 0); | |
| 5164 | 656 } |
| 657 | |
| 658 SparseComplexMatrix | |
| 5785 | 659 SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info) const |
| 5164 | 660 { |
| 661 double rcond; | |
| 5506 | 662 return inverse (mattype, info, rcond, 0, 0); |
| 663 } | |
| 664 | |
| 665 SparseComplexMatrix | |
| 5785 | 666 SparseComplexMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info, |
| 5610 | 667 double& rcond, const bool, |
| 5506 | 668 const bool calccond) const |
| 669 { | |
| 670 SparseComplexMatrix retval; | |
| 671 | |
| 672 octave_idx_type nr = rows (); | |
| 673 octave_idx_type nc = cols (); | |
| 674 info = 0; | |
| 675 | |
| 676 if (nr == 0 || nc == 0 || nr != nc) | |
| 677 (*current_liboctave_error_handler) ("inverse requires square matrix"); | |
| 678 else | |
| 679 { | |
| 680 // Print spparms("spumoni") info if requested | |
| 681 int typ = mattyp.type (); | |
| 682 mattyp.info (); | |
| 683 | |
| 5785 | 684 if (typ == MatrixType::Diagonal || |
| 685 typ == MatrixType::Permuted_Diagonal) | |
| 5506 | 686 { |
| 5785 | 687 if (typ == MatrixType::Permuted_Diagonal) |
| 5506 | 688 retval = transpose(); |
| 689 else | |
| 690 retval = *this; | |
| 691 | |
| 692 // Force make_unique to be called | |
| 693 Complex *v = retval.data(); | |
| 694 | |
| 695 if (calccond) | |
| 696 { | |
| 697 double dmax = 0., dmin = octave_Inf; | |
| 698 for (octave_idx_type i = 0; i < nr; i++) | |
| 699 { | |
| 700 double tmp = std::abs(v[i]); | |
| 701 if (tmp > dmax) | |
| 702 dmax = tmp; | |
| 703 if (tmp < dmin) | |
| 704 dmin = tmp; | |
| 705 } | |
| 706 rcond = dmin / dmax; | |
| 707 } | |
| 708 | |
| 709 for (octave_idx_type i = 0; i < nr; i++) | |
| 710 v[i] = 1.0 / v[i]; | |
| 711 } | |
| 712 else | |
| 713 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 714 } | |
| 715 | |
| 716 return retval; | |
| 717 } | |
| 718 | |
| 719 SparseComplexMatrix | |
| 5785 | 720 SparseComplexMatrix::tinverse (MatrixType &mattyp, octave_idx_type& info, |
| 5610 | 721 double& rcond, const bool, |
| 5506 | 722 const bool calccond) const |
| 723 { | |
| 724 SparseComplexMatrix retval; | |
| 725 | |
| 726 octave_idx_type nr = rows (); | |
| 727 octave_idx_type nc = cols (); | |
| 728 info = 0; | |
| 729 | |
| 730 if (nr == 0 || nc == 0 || nr != nc) | |
| 731 (*current_liboctave_error_handler) ("inverse requires square matrix"); | |
| 732 else | |
| 733 { | |
| 734 // Print spparms("spumoni") info if requested | |
| 735 int typ = mattyp.type (); | |
| 736 mattyp.info (); | |
| 737 | |
| 5785 | 738 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper || |
| 739 typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) | |
| 5506 | 740 { |
| 741 double anorm = 0.; | |
| 742 double ainvnorm = 0.; | |
| 743 | |
| 744 if (calccond) | |
| 745 { | |
| 746 // Calculate the 1-norm of matrix for rcond calculation | |
| 747 for (octave_idx_type j = 0; j < nr; j++) | |
| 748 { | |
| 749 double atmp = 0.; | |
| 750 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 751 atmp += std::abs(data(i)); | |
| 752 if (atmp > anorm) | |
| 753 anorm = atmp; | |
| 754 } | |
| 755 } | |
| 756 | |
| 5785 | 757 if (typ == MatrixType::Upper || typ == MatrixType::Lower) |
| 5506 | 758 { |
| 5681 | 759 octave_idx_type nz = nnz (); |
| 5506 | 760 octave_idx_type cx = 0; |
| 761 octave_idx_type nz2 = nz; | |
| 762 retval = SparseComplexMatrix (nr, nc, nz2); | |
| 763 | |
| 764 for (octave_idx_type i = 0; i < nr; i++) | |
| 765 { | |
| 766 OCTAVE_QUIT; | |
| 767 // place the 1 in the identity position | |
| 768 octave_idx_type cx_colstart = cx; | |
| 769 | |
| 770 if (cx == nz2) | |
| 771 { | |
| 772 nz2 *= 2; | |
| 773 retval.change_capacity (nz2); | |
| 774 } | |
| 775 | |
| 776 retval.xcidx(i) = cx; | |
| 777 retval.xridx(cx) = i; | |
| 778 retval.xdata(cx) = 1.0; | |
| 779 cx++; | |
| 780 | |
| 781 // iterate accross columns of input matrix | |
| 782 for (octave_idx_type j = i+1; j < nr; j++) | |
| 783 { | |
| 784 Complex v = 0.; | |
| 785 // iterate to calculate sum | |
| 786 octave_idx_type colXp = retval.xcidx(i); | |
| 787 octave_idx_type colUp = cidx(j); | |
| 788 octave_idx_type rpX, rpU; | |
| 5876 | 789 |
| 790 if (cidx(j) == cidx(j+1)) | |
| 791 { | |
| 792 (*current_liboctave_error_handler) | |
| 793 ("division by zero"); | |
| 794 goto inverse_singular; | |
| 795 } | |
| 796 | |
| 5506 | 797 do |
| 798 { | |
| 799 OCTAVE_QUIT; | |
| 800 rpX = retval.xridx(colXp); | |
| 801 rpU = ridx(colUp); | |
| 802 | |
| 803 if (rpX < rpU) | |
| 804 colXp++; | |
| 805 else if (rpX > rpU) | |
| 806 colUp++; | |
| 807 else | |
| 808 { | |
| 809 v -= retval.xdata(colXp) * data(colUp); | |
| 810 colXp++; | |
| 811 colUp++; | |
| 812 } | |
| 813 } while ((rpX<j) && (rpU<j) && | |
| 814 (colXp<cx) && (colUp<nz)); | |
| 815 | |
| 5876 | 816 |
| 5506 | 817 // get A(m,m) |
| 5876 | 818 if (typ == MatrixType::Upper) |
| 819 colUp = cidx(j+1) - 1; | |
| 820 else | |
| 5877 | 821 colUp = cidx(j); |
| 5506 | 822 Complex pivot = data(colUp); |
| 5877 | 823 if (pivot == 0. || ridx(colUp) != j) |
| 5876 | 824 { |
| 825 (*current_liboctave_error_handler) | |
| 826 ("division by zero"); | |
| 827 goto inverse_singular; | |
| 828 } | |
| 5506 | 829 |
| 830 if (v != 0.) | |
| 831 { | |
| 832 if (cx == nz2) | |
| 833 { | |
| 834 nz2 *= 2; | |
| 835 retval.change_capacity (nz2); | |
| 836 } | |
| 837 | |
| 838 retval.xridx(cx) = j; | |
| 839 retval.xdata(cx) = v / pivot; | |
| 840 cx++; | |
| 841 } | |
| 842 } | |
| 843 | |
| 844 // get A(m,m) | |
| 5876 | 845 octave_idx_type colUp; |
| 846 if (typ == MatrixType::Upper) | |
| 847 colUp = cidx(i+1) - 1; | |
| 848 else | |
| 5877 | 849 colUp = cidx(i); |
| 5506 | 850 Complex pivot = data(colUp); |
| 5877 | 851 if (pivot == 0. || ridx(colUp) != i) |
| 5876 | 852 { |
| 853 (*current_liboctave_error_handler) ("division by zero"); | |
| 854 goto inverse_singular; | |
| 855 } | |
| 5506 | 856 |
| 857 if (pivot != 1.0) | |
| 858 for (octave_idx_type j = cx_colstart; j < cx; j++) | |
| 859 retval.xdata(j) /= pivot; | |
| 860 } | |
| 861 retval.xcidx(nr) = cx; | |
| 862 retval.maybe_compress (); | |
| 863 } | |
| 864 else | |
| 865 { | |
| 5681 | 866 octave_idx_type nz = nnz (); |
| 5506 | 867 octave_idx_type cx = 0; |
| 868 octave_idx_type nz2 = nz; | |
| 869 retval = SparseComplexMatrix (nr, nc, nz2); | |
| 870 | |
| 871 OCTAVE_LOCAL_BUFFER (Complex, work, nr); | |
| 872 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); | |
| 873 | |
| 874 octave_idx_type *perm = mattyp.triangular_perm(); | |
| 5785 | 875 if (typ == MatrixType::Permuted_Upper) |
| 5506 | 876 { |
| 877 for (octave_idx_type i = 0; i < nr; i++) | |
| 878 rperm[perm[i]] = i; | |
| 879 } | |
| 880 else | |
| 881 { | |
| 882 for (octave_idx_type i = 0; i < nr; i++) | |
| 883 rperm[i] = perm[i]; | |
| 884 for (octave_idx_type i = 0; i < nr; i++) | |
| 885 perm[rperm[i]] = i; | |
| 886 } | |
| 887 | |
| 888 for (octave_idx_type i = 0; i < nr; i++) | |
| 889 { | |
| 890 OCTAVE_QUIT; | |
| 891 octave_idx_type iidx = rperm[i]; | |
| 892 | |
| 893 for (octave_idx_type j = 0; j < nr; j++) | |
| 894 work[j] = 0.; | |
| 895 | |
| 896 // place the 1 in the identity position | |
| 897 work[iidx] = 1.0; | |
| 898 | |
| 899 // iterate accross columns of input matrix | |
| 900 for (octave_idx_type j = iidx+1; j < nr; j++) | |
| 901 { | |
| 902 Complex v = 0.; | |
| 903 octave_idx_type jidx = perm[j]; | |
| 904 // iterate to calculate sum | |
| 905 for (octave_idx_type k = cidx(jidx); | |
| 906 k < cidx(jidx+1); k++) | |
| 907 { | |
| 908 OCTAVE_QUIT; | |
| 909 v -= work[ridx(k)] * data(k); | |
| 910 } | |
| 911 | |
| 912 // get A(m,m) | |
| 5876 | 913 Complex pivot; |
| 914 if (typ == MatrixType::Permuted_Upper) | |
| 915 pivot = data(cidx(jidx+1) - 1); | |
| 916 else | |
| 5877 | 917 pivot = data(cidx(jidx)); |
| 5506 | 918 if (pivot == 0.) |
| 5876 | 919 { |
| 920 (*current_liboctave_error_handler) | |
| 921 ("division by zero"); | |
| 922 goto inverse_singular; | |
| 923 } | |
| 5506 | 924 |
| 925 work[j] = v / pivot; | |
| 926 } | |
| 927 | |
| 928 // get A(m,m) | |
| 5876 | 929 octave_idx_type colUp; |
| 930 if (typ == MatrixType::Permuted_Upper) | |
| 931 colUp = cidx(perm[iidx]+1) - 1; | |
| 932 else | |
| 5877 | 933 colUp = cidx(perm[iidx]); |
| 5876 | 934 |
| 935 Complex pivot = data(colUp); | |
| 936 if (pivot == 0.) | |
| 937 { | |
| 938 (*current_liboctave_error_handler) | |
| 939 ("division by zero"); | |
| 940 goto inverse_singular; | |
| 941 } | |
| 5506 | 942 |
| 943 octave_idx_type new_cx = cx; | |
| 944 for (octave_idx_type j = iidx; j < nr; j++) | |
| 945 if (work[j] != 0.0) | |
| 946 { | |
| 947 new_cx++; | |
| 948 if (pivot != 1.0) | |
| 949 work[j] /= pivot; | |
| 950 } | |
| 951 | |
| 952 if (cx < new_cx) | |
| 953 { | |
| 954 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2); | |
| 955 retval.change_capacity (nz2); | |
| 956 } | |
| 957 | |
| 958 retval.xcidx(i) = cx; | |
| 959 for (octave_idx_type j = iidx; j < nr; j++) | |
| 960 if (work[j] != 0.) | |
| 961 { | |
| 962 retval.xridx(cx) = j; | |
| 963 retval.xdata(cx++) = work[j]; | |
| 964 } | |
| 965 } | |
| 966 | |
| 967 retval.xcidx(nr) = cx; | |
| 968 retval.maybe_compress (); | |
| 969 } | |
| 970 | |
| 971 if (calccond) | |
| 972 { | |
| 973 // Calculate the 1-norm of inverse matrix for rcond calculation | |
| 974 for (octave_idx_type j = 0; j < nr; j++) | |
| 975 { | |
| 976 double atmp = 0.; | |
| 977 for (octave_idx_type i = retval.cidx(j); | |
| 978 i < retval.cidx(j+1); i++) | |
| 979 atmp += std::abs(retval.data(i)); | |
| 980 if (atmp > ainvnorm) | |
| 981 ainvnorm = atmp; | |
| 982 } | |
| 983 | |
| 984 rcond = 1. / ainvnorm / anorm; | |
| 985 } | |
| 986 } | |
| 987 else | |
| 988 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 989 } | |
| 990 | |
| 991 return retval; | |
| 5876 | 992 |
| 993 inverse_singular: | |
| 994 return SparseComplexMatrix(); | |
| 5164 | 995 } |
| 996 | |
| 997 SparseComplexMatrix | |
| 5785 | 998 SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info, |
| 5610 | 999 double& rcond, int, int calc_cond) const |
| 5506 | 1000 { |
| 1001 int typ = mattype.type (false); | |
| 1002 SparseComplexMatrix ret; | |
| 1003 | |
| 5785 | 1004 if (typ == MatrixType::Unknown) |
| 5506 | 1005 typ = mattype.type (*this); |
| 1006 | |
| 5785 | 1007 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) |
| 5506 | 1008 ret = dinverse (mattype, info, rcond, true, calc_cond); |
| 5785 | 1009 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) |
| 5506 | 1010 ret = tinverse (mattype, info, rcond, true, calc_cond).transpose(); |
| 5785 | 1011 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) |
| 6185 | 1012 { |
| 1013 MatrixType newtype = mattype.transpose(); | |
| 1014 ret = transpose().tinverse (newtype, info, rcond, true, calc_cond); | |
| 1015 } | |
| 6840 | 1016 else |
| 5506 | 1017 { |
| 1018 if (mattype.is_hermitian()) | |
| 1019 { | |
| 5785 | 1020 MatrixType tmp_typ (MatrixType::Upper); |
| 5506 | 1021 SparseComplexCHOL fact (*this, info, false); |
| 1022 rcond = fact.rcond(); | |
| 1023 if (info == 0) | |
| 1024 { | |
| 1025 double rcond2; | |
| 1026 SparseMatrix Q = fact.Q(); | |
| 1027 SparseComplexMatrix InvL = fact.L().transpose(). | |
| 1028 tinverse(tmp_typ, info, rcond2, true, false); | |
| 1029 ret = Q * InvL.hermitian() * InvL * Q.transpose(); | |
| 1030 } | |
| 1031 else | |
| 1032 { | |
| 1033 // Matrix is either singular or not positive definite | |
| 1034 mattype.mark_as_unsymmetric (); | |
| 5785 | 1035 typ = MatrixType::Full; |
| 5506 | 1036 } |
| 1037 } | |
| 1038 | |
| 1039 if (!mattype.is_hermitian()) | |
| 1040 { | |
| 1041 octave_idx_type n = rows(); | |
| 1042 ColumnVector Qinit(n); | |
| 1043 for (octave_idx_type i = 0; i < n; i++) | |
| 1044 Qinit(i) = i; | |
| 1045 | |
| 5785 | 1046 MatrixType tmp_typ (MatrixType::Upper); |
| 5506 | 1047 SparseComplexLU fact (*this, Qinit, -1.0, false); |
| 1048 rcond = fact.rcond(); | |
| 1049 double rcond2; | |
| 1050 SparseComplexMatrix InvL = fact.L().transpose(). | |
| 1051 tinverse(tmp_typ, info, rcond2, true, false); | |
| 1052 SparseComplexMatrix InvU = fact.U(). | |
| 1053 tinverse(tmp_typ, info, rcond2, true, false).transpose(); | |
| 1054 ret = fact.Pc().transpose() * InvU * InvL * fact.Pr(); | |
| 1055 } | |
| 1056 } | |
| 1057 | |
| 1058 return ret; | |
| 5164 | 1059 } |
| 1060 | |
| 1061 ComplexDET | |
| 1062 SparseComplexMatrix::determinant (void) const | |
| 1063 { | |
| 5275 | 1064 octave_idx_type info; |
| 5164 | 1065 double rcond; |
| 1066 return determinant (info, rcond, 0); | |
| 1067 } | |
| 1068 | |
| 1069 ComplexDET | |
| 5275 | 1070 SparseComplexMatrix::determinant (octave_idx_type& info) const |
| 5164 | 1071 { |
| 1072 double rcond; | |
| 1073 return determinant (info, rcond, 0); | |
| 1074 } | |
| 1075 | |
| 1076 ComplexDET | |
| 5610 | 1077 SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, int) const |
| 5164 | 1078 { |
| 1079 ComplexDET retval; | |
| 5203 | 1080 #ifdef HAVE_UMFPACK |
| 5164 | 1081 |
| 5275 | 1082 octave_idx_type nr = rows (); |
| 1083 octave_idx_type nc = cols (); | |
| 5164 | 1084 |
| 1085 if (nr == 0 || nc == 0 || nr != nc) | |
| 1086 { | |
| 1087 Complex d[2]; | |
| 1088 d[0] = 1.0; | |
| 1089 d[1] = 0.0; | |
| 1090 retval = ComplexDET (d); | |
| 1091 } | |
| 1092 else | |
| 1093 { | |
| 1094 err = 0; | |
| 1095 | |
| 1096 // Setup the control parameters | |
| 1097 Matrix Control (UMFPACK_CONTROL, 1); | |
| 1098 double *control = Control.fortran_vec (); | |
| 5322 | 1099 UMFPACK_ZNAME (defaults) (control); |
| 5164 | 1100 |
| 5893 | 1101 double tmp = octave_sparse_params::get_key ("spumoni"); |
| 5164 | 1102 if (!xisnan (tmp)) |
| 1103 Control (UMFPACK_PRL) = tmp; | |
| 1104 | |
| 5893 | 1105 tmp = octave_sparse_params::get_key ("piv_tol"); |
| 5164 | 1106 if (!xisnan (tmp)) |
| 1107 { | |
| 1108 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; | |
| 1109 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; | |
| 1110 } | |
| 1111 | |
| 1112 // Set whether we are allowed to modify Q or not | |
| 5893 | 1113 tmp = octave_sparse_params::get_key ("autoamd"); |
| 5164 | 1114 if (!xisnan (tmp)) |
| 1115 Control (UMFPACK_FIXQ) = tmp; | |
| 1116 | |
| 1117 // Turn-off UMFPACK scaling for LU | |
| 1118 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; | |
| 1119 | |
| 5322 | 1120 UMFPACK_ZNAME (report_control) (control); |
| 5164 | 1121 |
| 5275 | 1122 const octave_idx_type *Ap = cidx (); |
| 1123 const octave_idx_type *Ai = ridx (); | |
| 5164 | 1124 const Complex *Ax = data (); |
| 1125 | |
| 5322 | 1126 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, |
| 5760 | 1127 reinterpret_cast<const double *> (Ax), |
| 1128 NULL, 1, control); | |
| 5164 | 1129 |
| 1130 void *Symbolic; | |
| 1131 Matrix Info (1, UMFPACK_INFO); | |
| 1132 double *info = Info.fortran_vec (); | |
| 5322 | 1133 int status = UMFPACK_ZNAME (qsymbolic) |
| 5760 | 1134 (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), NULL, |
| 5164 | 1135 NULL, &Symbolic, control, info); |
| 1136 | |
| 1137 if (status < 0) | |
| 1138 { | |
| 1139 (*current_liboctave_error_handler) | |
| 1140 ("SparseComplexMatrix::determinant symbolic factorization failed"); | |
| 1141 | |
| 5322 | 1142 UMFPACK_ZNAME (report_status) (control, status); |
| 1143 UMFPACK_ZNAME (report_info) (control, info); | |
| 1144 | |
| 1145 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; | |
| 5164 | 1146 } |
| 1147 else | |
| 1148 { | |
| 5322 | 1149 UMFPACK_ZNAME (report_symbolic) (Symbolic, control); |
| 5164 | 1150 |
| 1151 void *Numeric; | |
| 5760 | 1152 status |
| 1153 = UMFPACK_ZNAME (numeric) (Ap, Ai, | |
| 1154 reinterpret_cast<const double *> (Ax), | |
| 1155 NULL, Symbolic, &Numeric, control, info) ; | |
| 5322 | 1156 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; |
| 5164 | 1157 |
| 1158 rcond = Info (UMFPACK_RCOND); | |
| 1159 | |
| 1160 if (status < 0) | |
| 1161 { | |
| 1162 (*current_liboctave_error_handler) | |
| 1163 ("SparseComplexMatrix::determinant numeric factorization failed"); | |
| 1164 | |
| 5322 | 1165 UMFPACK_ZNAME (report_status) (control, status); |
| 1166 UMFPACK_ZNAME (report_info) (control, info); | |
| 1167 | |
| 1168 UMFPACK_ZNAME (free_numeric) (&Numeric); | |
| 5164 | 1169 } |
| 1170 else | |
| 1171 { | |
| 5322 | 1172 UMFPACK_ZNAME (report_numeric) (Numeric, control); |
| 5164 | 1173 |
| 1174 Complex d[2]; | |
| 1175 double d_exponent; | |
| 1176 | |
| 5322 | 1177 status = UMFPACK_ZNAME (get_determinant) |
| 5760 | 1178 (reinterpret_cast<double *> (&d[0]), NULL, &d_exponent, |
| 5164 | 1179 Numeric, info); |
| 1180 d[1] = d_exponent; | |
| 1181 | |
| 1182 if (status < 0) | |
| 1183 { | |
| 1184 (*current_liboctave_error_handler) | |
| 1185 ("SparseComplexMatrix::determinant error calculating determinant"); | |
| 1186 | |
| 5322 | 1187 UMFPACK_ZNAME (report_status) (control, status); |
| 1188 UMFPACK_ZNAME (report_info) (control, info); | |
| 5164 | 1189 } |
| 1190 else | |
| 1191 retval = ComplexDET (d); | |
| 5346 | 1192 |
| 1193 UMFPACK_ZNAME (free_numeric) (&Numeric); | |
| 5164 | 1194 } |
| 1195 } | |
| 1196 } | |
| 5203 | 1197 #else |
| 1198 (*current_liboctave_error_handler) ("UMFPACK not installed"); | |
| 1199 #endif | |
| 5164 | 1200 |
| 1201 return retval; | |
| 1202 } | |
| 1203 | |
| 1204 ComplexMatrix | |
| 5785 | 1205 SparseComplexMatrix::dsolve (MatrixType &mattype, const Matrix& b, |
| 5681 | 1206 octave_idx_type& err, double& rcond, |
| 1207 solve_singularity_handler, bool calc_cond) const | |
| 5164 | 1208 { |
| 1209 ComplexMatrix retval; | |
| 1210 | |
| 5275 | 1211 octave_idx_type nr = rows (); |
| 1212 octave_idx_type nc = cols (); | |
| 5630 | 1213 octave_idx_type nm = (nc < nr ? nc : nr); |
| 5164 | 1214 err = 0; |
| 1215 | |
| 6924 | 1216 if (nr != b.rows ()) |
| 5164 | 1217 (*current_liboctave_error_handler) |
| 1218 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 1219 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 1220 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
| 5164 | 1221 else |
| 1222 { | |
| 1223 // Print spparms("spumoni") info if requested | |
| 1224 int typ = mattype.type (); | |
| 1225 mattype.info (); | |
| 1226 | |
| 5785 | 1227 if (typ == MatrixType::Diagonal || |
| 1228 typ == MatrixType::Permuted_Diagonal) | |
| 5164 | 1229 { |
| 5630 | 1230 retval.resize (nc, b.cols(), Complex(0.,0.)); |
| 5785 | 1231 if (typ == MatrixType::Diagonal) |
| 5275 | 1232 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5630 | 1233 for (octave_idx_type i = 0; i < nm; i++) |
| 1234 retval(i,j) = b(i,j) / data (i); | |
| 5164 | 1235 else |
| 5275 | 1236 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5630 | 1237 for (octave_idx_type k = 0; k < nc; k++) |
| 1238 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) | |
| 1239 retval(k,j) = b(ridx(i),j) / data (i); | |
| 5164 | 1240 |
| 5681 | 1241 if (calc_cond) |
| 1242 { | |
| 1243 double dmax = 0., dmin = octave_Inf; | |
| 1244 for (octave_idx_type i = 0; i < nm; i++) | |
| 1245 { | |
| 1246 double tmp = std::abs(data(i)); | |
| 1247 if (tmp > dmax) | |
| 1248 dmax = tmp; | |
| 1249 if (tmp < dmin) | |
| 1250 dmin = tmp; | |
| 1251 } | |
| 1252 rcond = dmin / dmax; | |
| 1253 } | |
| 1254 else | |
| 1255 rcond = 1.0; | |
| 5164 | 1256 } |
| 1257 else | |
| 1258 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 1259 } | |
| 1260 | |
| 1261 return retval; | |
| 1262 } | |
| 1263 | |
| 1264 SparseComplexMatrix | |
| 5785 | 1265 SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseMatrix& b, |
| 5630 | 1266 octave_idx_type& err, double& rcond, |
| 5681 | 1267 solve_singularity_handler, |
| 1268 bool calc_cond) const | |
| 5164 | 1269 { |
| 1270 SparseComplexMatrix retval; | |
| 1271 | |
| 5275 | 1272 octave_idx_type nr = rows (); |
| 1273 octave_idx_type nc = cols (); | |
| 5630 | 1274 octave_idx_type nm = (nc < nr ? nc : nr); |
| 5164 | 1275 err = 0; |
| 1276 | |
| 6924 | 1277 if (nr != b.rows ()) |
| 5164 | 1278 (*current_liboctave_error_handler) |
| 1279 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 1280 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 1281 retval = SparseComplexMatrix (nc, b.cols ()); | |
| 5164 | 1282 else |
| 1283 { | |
| 1284 // Print spparms("spumoni") info if requested | |
| 1285 int typ = mattype.type (); | |
| 1286 mattype.info (); | |
| 1287 | |
| 5785 | 1288 if (typ == MatrixType::Diagonal || |
| 1289 typ == MatrixType::Permuted_Diagonal) | |
| 5164 | 1290 { |
| 5275 | 1291 octave_idx_type b_nc = b.cols (); |
| 5681 | 1292 octave_idx_type b_nz = b.nnz (); |
| 5630 | 1293 retval = SparseComplexMatrix (nc, b_nc, b_nz); |
| 5164 | 1294 |
| 1295 retval.xcidx(0) = 0; | |
| 5275 | 1296 octave_idx_type ii = 0; |
| 5785 | 1297 if (typ == MatrixType::Diagonal) |
| 5275 | 1298 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5164 | 1299 { |
| 5275 | 1300 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) |
| 5164 | 1301 { |
| 5681 | 1302 if (b.ridx(i) >= nm) |
| 1303 break; | |
| 5164 | 1304 retval.xridx (ii) = b.ridx(i); |
| 1305 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); | |
| 1306 } | |
| 1307 retval.xcidx(j+1) = ii; | |
| 1308 } | |
| 1309 else | |
| 5275 | 1310 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5164 | 1311 { |
| 5630 | 1312 for (octave_idx_type l = 0; l < nc; l++) |
| 1313 for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) | |
| 1314 { | |
| 1315 bool found = false; | |
| 1316 octave_idx_type k; | |
| 1317 for (k = b.cidx(j); k < b.cidx(j+1); k++) | |
| 1318 if (ridx(i) == b.ridx(k)) | |
| 1319 { | |
| 1320 found = true; | |
| 1321 break; | |
| 1322 } | |
| 1323 if (found) | |
| 5164 | 1324 { |
| 5630 | 1325 retval.xridx (ii) = l; |
| 1326 retval.xdata (ii++) = b.data(k) / data (i); | |
| 5164 | 1327 } |
| 5630 | 1328 } |
| 5164 | 1329 retval.xcidx(j+1) = ii; |
| 1330 } | |
| 1331 | |
| 5681 | 1332 if (calc_cond) |
| 1333 { | |
| 1334 double dmax = 0., dmin = octave_Inf; | |
| 1335 for (octave_idx_type i = 0; i < nm; i++) | |
| 1336 { | |
| 1337 double tmp = std::abs(data(i)); | |
| 1338 if (tmp > dmax) | |
| 1339 dmax = tmp; | |
| 1340 if (tmp < dmin) | |
| 1341 dmin = tmp; | |
| 1342 } | |
| 1343 rcond = dmin / dmax; | |
| 1344 } | |
| 1345 else | |
| 1346 rcond = 1.0; | |
| 5164 | 1347 } |
| 1348 else | |
| 1349 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 1350 } | |
| 1351 | |
| 1352 return retval; | |
| 1353 } | |
| 1354 | |
| 1355 ComplexMatrix | |
| 5785 | 1356 SparseComplexMatrix::dsolve (MatrixType &mattype, const ComplexMatrix& b, |
| 5630 | 1357 octave_idx_type& err, double& rcond, |
| 5681 | 1358 solve_singularity_handler, |
| 1359 bool calc_cond) const | |
| 5164 | 1360 { |
| 1361 ComplexMatrix retval; | |
| 1362 | |
| 5275 | 1363 octave_idx_type nr = rows (); |
| 1364 octave_idx_type nc = cols (); | |
| 5630 | 1365 octave_idx_type nm = (nc < nr ? nc : nr); |
| 5164 | 1366 err = 0; |
| 1367 | |
| 6924 | 1368 if (nr != b.rows ()) |
| 5164 | 1369 (*current_liboctave_error_handler) |
| 1370 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 1371 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 1372 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
| 5164 | 1373 else |
| 1374 { | |
| 1375 // Print spparms("spumoni") info if requested | |
| 1376 int typ = mattype.type (); | |
| 1377 mattype.info (); | |
| 1378 | |
| 5785 | 1379 if (typ == MatrixType::Diagonal || |
| 1380 typ == MatrixType::Permuted_Diagonal) | |
| 5164 | 1381 { |
| 5630 | 1382 retval.resize (nc, b.cols(), Complex(0.,0.)); |
| 5785 | 1383 if (typ == MatrixType::Diagonal) |
| 5275 | 1384 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5630 | 1385 for (octave_idx_type i = 0; i < nm; i++) |
| 5164 | 1386 retval(i,j) = b(i,j) / data (i); |
| 1387 else | |
| 5275 | 1388 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5630 | 1389 for (octave_idx_type k = 0; k < nc; k++) |
| 1390 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) | |
| 1391 retval(k,j) = b(ridx(i),j) / data (i); | |
| 5164 | 1392 |
| 5681 | 1393 if (calc_cond) |
| 1394 { | |
| 1395 double dmax = 0., dmin = octave_Inf; | |
| 1396 for (octave_idx_type i = 0; i < nr; i++) | |
| 1397 { | |
| 1398 double tmp = std::abs(data(i)); | |
| 1399 if (tmp > dmax) | |
| 1400 dmax = tmp; | |
| 1401 if (tmp < dmin) | |
| 1402 dmin = tmp; | |
| 1403 } | |
| 1404 rcond = dmin / dmax; | |
| 1405 } | |
| 1406 else | |
| 1407 rcond = 1.0; | |
| 5164 | 1408 } |
| 1409 else | |
| 1410 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 1411 } | |
| 1412 | |
| 1413 return retval; | |
| 1414 } | |
| 1415 | |
| 1416 SparseComplexMatrix | |
| 5785 | 1417 SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseComplexMatrix& b, |
| 5630 | 1418 octave_idx_type& err, double& rcond, |
| 5681 | 1419 solve_singularity_handler, |
| 1420 bool calc_cond) const | |
| 5164 | 1421 { |
| 1422 SparseComplexMatrix retval; | |
| 1423 | |
| 5275 | 1424 octave_idx_type nr = rows (); |
| 1425 octave_idx_type nc = cols (); | |
| 5630 | 1426 octave_idx_type nm = (nc < nr ? nc : nr); |
| 5164 | 1427 err = 0; |
| 1428 | |
| 6924 | 1429 if (nr != b.rows ()) |
| 5164 | 1430 (*current_liboctave_error_handler) |
| 1431 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 1432 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 1433 retval = SparseComplexMatrix (nc, b.cols ()); | |
| 5164 | 1434 else |
| 1435 { | |
| 1436 // Print spparms("spumoni") info if requested | |
| 1437 int typ = mattype.type (); | |
| 1438 mattype.info (); | |
| 1439 | |
| 5785 | 1440 if (typ == MatrixType::Diagonal || |
| 1441 typ == MatrixType::Permuted_Diagonal) | |
| 5164 | 1442 { |
| 5275 | 1443 octave_idx_type b_nc = b.cols (); |
| 5681 | 1444 octave_idx_type b_nz = b.nnz (); |
| 5630 | 1445 retval = SparseComplexMatrix (nc, b_nc, b_nz); |
| 5164 | 1446 |
| 1447 retval.xcidx(0) = 0; | |
| 5275 | 1448 octave_idx_type ii = 0; |
| 5785 | 1449 if (typ == MatrixType::Diagonal) |
| 5275 | 1450 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5164 | 1451 { |
| 5275 | 1452 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) |
| 5164 | 1453 { |
| 5681 | 1454 if (b.ridx(i) >= nm) |
| 1455 break; | |
| 5164 | 1456 retval.xridx (ii) = b.ridx(i); |
| 1457 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); | |
| 1458 } | |
| 1459 retval.xcidx(j+1) = ii; | |
| 1460 } | |
| 1461 else | |
| 5275 | 1462 for (octave_idx_type j = 0; j < b.cols(); j++) |
| 5164 | 1463 { |
| 5630 | 1464 for (octave_idx_type l = 0; l < nc; l++) |
| 1465 for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) | |
| 1466 { | |
| 1467 bool found = false; | |
| 1468 octave_idx_type k; | |
| 1469 for (k = b.cidx(j); k < b.cidx(j+1); k++) | |
| 1470 if (ridx(i) == b.ridx(k)) | |
| 1471 { | |
| 1472 found = true; | |
| 1473 break; | |
| 1474 } | |
| 1475 if (found) | |
| 5164 | 1476 { |
| 5630 | 1477 retval.xridx (ii) = l; |
| 1478 retval.xdata (ii++) = b.data(k) / data (i); | |
| 5164 | 1479 } |
| 5630 | 1480 } |
| 5164 | 1481 retval.xcidx(j+1) = ii; |
| 1482 } | |
| 1483 | |
| 5681 | 1484 if (calc_cond) |
| 1485 { | |
| 1486 double dmax = 0., dmin = octave_Inf; | |
| 1487 for (octave_idx_type i = 0; i < nm; i++) | |
| 1488 { | |
| 1489 double tmp = std::abs(data(i)); | |
| 1490 if (tmp > dmax) | |
| 1491 dmax = tmp; | |
| 1492 if (tmp < dmin) | |
| 1493 dmin = tmp; | |
| 1494 } | |
| 1495 rcond = dmin / dmax; | |
| 1496 } | |
| 1497 else | |
| 1498 rcond = 1.0; | |
| 5164 | 1499 } |
| 1500 else | |
| 1501 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 1502 } | |
| 1503 | |
| 1504 return retval; | |
| 1505 } | |
| 1506 | |
| 1507 ComplexMatrix | |
| 5785 | 1508 SparseComplexMatrix::utsolve (MatrixType &mattype, const Matrix& b, |
| 5630 | 1509 octave_idx_type& err, double& rcond, |
| 5681 | 1510 solve_singularity_handler sing_handler, |
| 1511 bool calc_cond) const | |
| 5164 | 1512 { |
| 1513 ComplexMatrix retval; | |
| 1514 | |
| 5275 | 1515 octave_idx_type nr = rows (); |
| 1516 octave_idx_type nc = cols (); | |
| 5630 | 1517 octave_idx_type nm = (nc > nr ? nc : nr); |
| 5164 | 1518 err = 0; |
| 1519 | |
| 6924 | 1520 if (nr != b.rows ()) |
| 5164 | 1521 (*current_liboctave_error_handler) |
| 1522 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 1523 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 1524 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
| 5164 | 1525 else |
| 1526 { | |
| 1527 // Print spparms("spumoni") info if requested | |
| 1528 int typ = mattype.type (); | |
| 1529 mattype.info (); | |
| 1530 | |
| 5785 | 1531 if (typ == MatrixType::Permuted_Upper || |
| 1532 typ == MatrixType::Upper) | |
| 5164 | 1533 { |
| 1534 double anorm = 0.; | |
| 1535 double ainvnorm = 0.; | |
| 5630 | 1536 octave_idx_type b_nc = b.cols (); |
| 5681 | 1537 rcond = 1.; |
| 1538 | |
| 1539 if (calc_cond) | |
| 1540 { | |
| 1541 // Calculate the 1-norm of matrix for rcond calculation | |
| 1542 for (octave_idx_type j = 0; j < nc; j++) | |
| 1543 { | |
| 1544 double atmp = 0.; | |
| 1545 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 1546 atmp += std::abs(data(i)); | |
| 1547 if (atmp > anorm) | |
| 1548 anorm = atmp; | |
| 1549 } | |
| 5164 | 1550 } |
| 1551 | |
| 5785 | 1552 if (typ == MatrixType::Permuted_Upper) |
| 5164 | 1553 { |
| 5630 | 1554 retval.resize (nc, b_nc); |
| 5322 | 1555 octave_idx_type *perm = mattype.triangular_perm (); |
| 5681 | 1556 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 5164 | 1557 |
| 5630 | 1558 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 1559 { |
| 5275 | 1560 for (octave_idx_type i = 0; i < nr; i++) |
| 5164 | 1561 work[i] = b(i,j); |
| 5630 | 1562 for (octave_idx_type i = nr; i < nc; i++) |
| 1563 work[i] = 0.; | |
| 1564 | |
| 1565 for (octave_idx_type k = nc-1; k >= 0; k--) | |
| 5164 | 1566 { |
| 5322 | 1567 octave_idx_type kidx = perm[k]; |
| 1568 | |
| 1569 if (work[k] != 0.) | |
| 5164 | 1570 { |
| 5681 | 1571 if (ridx(cidx(kidx+1)-1) != k || |
| 1572 data(cidx(kidx+1)-1) == 0.) | |
| 5164 | 1573 { |
| 1574 err = -2; | |
| 1575 goto triangular_error; | |
| 1576 } | |
| 1577 | |
| 5322 | 1578 Complex tmp = work[k] / data(cidx(kidx+1)-1); |
| 1579 work[k] = tmp; | |
| 1580 for (octave_idx_type i = cidx(kidx); | |
| 1581 i < cidx(kidx+1)-1; i++) | |
| 5164 | 1582 { |
| 5322 | 1583 octave_idx_type iidx = ridx(i); |
| 1584 work[iidx] = work[iidx] - tmp * data(i); | |
| 5164 | 1585 } |
| 1586 } | |
| 1587 } | |
| 1588 | |
| 5630 | 1589 for (octave_idx_type i = 0; i < nc; i++) |
| 5322 | 1590 retval (perm[i], j) = work[i]; |
| 5164 | 1591 } |
| 1592 | |
| 5681 | 1593 if (calc_cond) |
| 1594 { | |
| 1595 // Calculation of 1-norm of inv(*this) | |
| 1596 for (octave_idx_type i = 0; i < nm; i++) | |
| 1597 work[i] = 0.; | |
| 1598 | |
| 1599 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 1600 { |
| 5681 | 1601 work[j] = 1.; |
| 1602 | |
| 1603 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 1604 { |
| 5681 | 1605 octave_idx_type iidx = perm[k]; |
| 1606 | |
| 1607 if (work[k] != 0.) | |
| 5164 | 1608 { |
| 5681 | 1609 Complex tmp = work[k] / data(cidx(iidx+1)-1); |
| 1610 work[k] = tmp; | |
| 1611 for (octave_idx_type i = cidx(iidx); | |
| 1612 i < cidx(iidx+1)-1; i++) | |
| 1613 { | |
| 1614 octave_idx_type idx2 = ridx(i); | |
| 1615 work[idx2] = work[idx2] - tmp * data(i); | |
| 1616 } | |
| 5164 | 1617 } |
| 1618 } | |
| 5681 | 1619 double atmp = 0; |
| 1620 for (octave_idx_type i = 0; i < j+1; i++) | |
| 1621 { | |
| 1622 atmp += std::abs(work[i]); | |
| 1623 work[i] = 0.; | |
| 1624 } | |
| 1625 if (atmp > ainvnorm) | |
| 1626 ainvnorm = atmp; | |
| 5164 | 1627 } |
| 5681 | 1628 rcond = 1. / ainvnorm / anorm; |
| 5164 | 1629 } |
| 1630 } | |
| 1631 else | |
| 1632 { | |
| 5630 | 1633 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 1634 retval.resize (nc, b_nc); | |
| 1635 | |
| 1636 for (octave_idx_type j = 0; j < b_nc; j++) | |
| 5164 | 1637 { |
| 5630 | 1638 for (octave_idx_type i = 0; i < nr; i++) |
| 1639 work[i] = b(i,j); | |
| 1640 for (octave_idx_type i = nr; i < nc; i++) | |
| 1641 work[i] = 0.; | |
| 1642 | |
| 1643 for (octave_idx_type k = nc-1; k >= 0; k--) | |
| 5164 | 1644 { |
| 5630 | 1645 if (work[k] != 0.) |
| 5164 | 1646 { |
| 5681 | 1647 if (ridx(cidx(k+1)-1) != k || |
| 1648 data(cidx(k+1)-1) == 0.) | |
| 5164 | 1649 { |
| 1650 err = -2; | |
| 1651 goto triangular_error; | |
| 1652 } | |
| 1653 | |
| 5630 | 1654 Complex tmp = work[k] / data(cidx(k+1)-1); |
| 1655 work[k] = tmp; | |
| 5275 | 1656 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) |
| 5164 | 1657 { |
| 5275 | 1658 octave_idx_type iidx = ridx(i); |
| 5630 | 1659 work[iidx] = work[iidx] - tmp * data(i); |
| 5164 | 1660 } |
| 1661 } | |
| 1662 } | |
| 5630 | 1663 |
| 1664 for (octave_idx_type i = 0; i < nc; i++) | |
| 1665 retval.xelem (i, j) = work[i]; | |
| 5164 | 1666 } |
| 1667 | |
| 5681 | 1668 if (calc_cond) |
| 1669 { | |
| 1670 // Calculation of 1-norm of inv(*this) | |
| 1671 for (octave_idx_type i = 0; i < nm; i++) | |
| 1672 work[i] = 0.; | |
| 1673 | |
| 1674 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 1675 { |
| 5681 | 1676 work[j] = 1.; |
| 1677 | |
| 1678 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 1679 { |
| 5681 | 1680 if (work[k] != 0.) |
| 5164 | 1681 { |
| 5681 | 1682 Complex tmp = work[k] / data(cidx(k+1)-1); |
| 1683 work[k] = tmp; | |
| 1684 for (octave_idx_type i = cidx(k); | |
| 1685 i < cidx(k+1)-1; i++) | |
| 1686 { | |
| 1687 octave_idx_type iidx = ridx(i); | |
| 1688 work[iidx] = work[iidx] - tmp * data(i); | |
| 1689 } | |
| 5164 | 1690 } |
| 1691 } | |
| 5681 | 1692 double atmp = 0; |
| 1693 for (octave_idx_type i = 0; i < j+1; i++) | |
| 1694 { | |
| 1695 atmp += std::abs(work[i]); | |
| 1696 work[i] = 0.; | |
| 1697 } | |
| 1698 if (atmp > ainvnorm) | |
| 1699 ainvnorm = atmp; | |
| 5164 | 1700 } |
| 5681 | 1701 rcond = 1. / ainvnorm / anorm; |
| 1702 } | |
| 1703 } | |
| 5164 | 1704 |
| 1705 triangular_error: | |
| 1706 if (err != 0) | |
| 1707 { | |
| 1708 if (sing_handler) | |
| 5681 | 1709 { |
| 1710 sing_handler (rcond); | |
| 1711 mattype.mark_as_rectangular (); | |
| 1712 } | |
| 5164 | 1713 else |
| 1714 (*current_liboctave_error_handler) | |
| 1715 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", | |
| 1716 rcond); | |
| 1717 } | |
| 1718 | |
| 1719 volatile double rcond_plus_one = rcond + 1.0; | |
| 1720 | |
| 1721 if (rcond_plus_one == 1.0 || xisnan (rcond)) | |
| 1722 { | |
| 1723 err = -2; | |
| 1724 | |
| 1725 if (sing_handler) | |
| 5681 | 1726 { |
| 1727 sing_handler (rcond); | |
| 1728 mattype.mark_as_rectangular (); | |
| 1729 } | |
| 5164 | 1730 else |
| 1731 (*current_liboctave_error_handler) | |
| 1732 ("matrix singular to machine precision, rcond = %g", | |
| 1733 rcond); | |
| 1734 } | |
| 1735 } | |
| 1736 else | |
| 1737 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 1738 } | |
| 1739 | |
| 1740 return retval; | |
| 1741 } | |
| 1742 | |
| 1743 SparseComplexMatrix | |
| 5785 | 1744 SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseMatrix& b, |
| 5630 | 1745 octave_idx_type& err, double& rcond, |
| 5681 | 1746 solve_singularity_handler sing_handler, |
| 1747 bool calc_cond) const | |
| 5164 | 1748 { |
| 1749 SparseComplexMatrix retval; | |
| 1750 | |
| 5275 | 1751 octave_idx_type nr = rows (); |
| 1752 octave_idx_type nc = cols (); | |
| 5630 | 1753 octave_idx_type nm = (nc > nr ? nc : nr); |
| 5164 | 1754 err = 0; |
| 1755 | |
| 6924 | 1756 if (nr != b.rows ()) |
| 5164 | 1757 (*current_liboctave_error_handler) |
| 1758 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 1759 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 1760 retval = SparseComplexMatrix (nc, b.cols ()); | |
| 5164 | 1761 else |
| 1762 { | |
| 1763 // Print spparms("spumoni") info if requested | |
| 1764 int typ = mattype.type (); | |
| 1765 mattype.info (); | |
| 1766 | |
| 5785 | 1767 if (typ == MatrixType::Permuted_Upper || |
| 1768 typ == MatrixType::Upper) | |
| 5164 | 1769 { |
| 1770 double anorm = 0.; | |
| 1771 double ainvnorm = 0.; | |
| 5681 | 1772 rcond = 1.; |
| 1773 | |
| 1774 if (calc_cond) | |
| 1775 { | |
| 1776 // Calculate the 1-norm of matrix for rcond calculation | |
| 1777 for (octave_idx_type j = 0; j < nc; j++) | |
| 1778 { | |
| 1779 double atmp = 0.; | |
| 1780 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 1781 atmp += std::abs(data(i)); | |
| 1782 if (atmp > anorm) | |
| 1783 anorm = atmp; | |
| 1784 } | |
| 5164 | 1785 } |
| 1786 | |
| 5275 | 1787 octave_idx_type b_nc = b.cols (); |
| 5681 | 1788 octave_idx_type b_nz = b.nnz (); |
| 5630 | 1789 retval = SparseComplexMatrix (nc, b_nc, b_nz); |
| 5164 | 1790 retval.xcidx(0) = 0; |
| 5275 | 1791 octave_idx_type ii = 0; |
| 1792 octave_idx_type x_nz = b_nz; | |
| 5164 | 1793 |
| 5785 | 1794 if (typ == MatrixType::Permuted_Upper) |
| 5164 | 1795 { |
| 5322 | 1796 octave_idx_type *perm = mattype.triangular_perm (); |
| 5630 | 1797 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 1798 | |
| 1799 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); | |
| 1800 for (octave_idx_type i = 0; i < nc; i++) | |
| 5322 | 1801 rperm[perm[i]] = i; |
| 5164 | 1802 |
| 5275 | 1803 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 1804 { |
| 5630 | 1805 for (octave_idx_type i = 0; i < nm; i++) |
| 5164 | 1806 work[i] = 0.; |
| 5275 | 1807 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) |
| 5164 | 1808 work[b.ridx(i)] = b.data(i); |
| 1809 | |
| 5630 | 1810 for (octave_idx_type k = nc-1; k >= 0; k--) |
| 5164 | 1811 { |
| 5322 | 1812 octave_idx_type kidx = perm[k]; |
| 1813 | |
| 1814 if (work[k] != 0.) | |
| 5164 | 1815 { |
| 5681 | 1816 if (ridx(cidx(kidx+1)-1) != k || |
| 1817 data(cidx(kidx+1)-1) == 0.) | |
| 5164 | 1818 { |
| 1819 err = -2; | |
| 1820 goto triangular_error; | |
| 1821 } | |
| 1822 | |
| 5322 | 1823 Complex tmp = work[k] / data(cidx(kidx+1)-1); |
| 1824 work[k] = tmp; | |
| 1825 for (octave_idx_type i = cidx(kidx); | |
| 1826 i < cidx(kidx+1)-1; i++) | |
| 5164 | 1827 { |
| 5322 | 1828 octave_idx_type iidx = ridx(i); |
| 1829 work[iidx] = work[iidx] - tmp * data(i); | |
| 5164 | 1830 } |
| 1831 } | |
| 1832 } | |
| 1833 | |
| 1834 // Count non-zeros in work vector and adjust space in | |
| 1835 // retval if needed | |
| 5275 | 1836 octave_idx_type new_nnz = 0; |
| 5630 | 1837 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 1838 if (work[i] != 0.) |
| 1839 new_nnz++; | |
| 1840 | |
| 1841 if (ii + new_nnz > x_nz) | |
| 1842 { | |
| 1843 // Resize the sparse matrix | |
| 5275 | 1844 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; |
| 5164 | 1845 retval.change_capacity (sz); |
| 1846 x_nz = sz; | |
| 1847 } | |
| 1848 | |
| 5630 | 1849 for (octave_idx_type i = 0; i < nc; i++) |
| 5322 | 1850 if (work[rperm[i]] != 0.) |
| 5164 | 1851 { |
| 1852 retval.xridx(ii) = i; | |
| 5322 | 1853 retval.xdata(ii++) = work[rperm[i]]; |
| 5164 | 1854 } |
| 1855 retval.xcidx(j+1) = ii; | |
| 1856 } | |
| 1857 | |
| 1858 retval.maybe_compress (); | |
| 1859 | |
| 5681 | 1860 if (calc_cond) |
| 1861 { | |
| 1862 // Calculation of 1-norm of inv(*this) | |
| 1863 for (octave_idx_type i = 0; i < nm; i++) | |
| 1864 work[i] = 0.; | |
| 1865 | |
| 1866 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 1867 { |
| 5681 | 1868 work[j] = 1.; |
| 1869 | |
| 1870 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 1871 { |
| 5681 | 1872 octave_idx_type iidx = perm[k]; |
| 1873 | |
| 1874 if (work[k] != 0.) | |
| 5164 | 1875 { |
| 5681 | 1876 Complex tmp = work[k] / data(cidx(iidx+1)-1); |
| 1877 work[k] = tmp; | |
| 1878 for (octave_idx_type i = cidx(iidx); | |
| 1879 i < cidx(iidx+1)-1; i++) | |
| 1880 { | |
| 1881 octave_idx_type idx2 = ridx(i); | |
| 1882 work[idx2] = work[idx2] - tmp * data(i); | |
| 1883 } | |
| 5164 | 1884 } |
| 1885 } | |
| 5681 | 1886 double atmp = 0; |
| 1887 for (octave_idx_type i = 0; i < j+1; i++) | |
| 1888 { | |
| 1889 atmp += std::abs(work[i]); | |
| 1890 work[i] = 0.; | |
| 1891 } | |
| 1892 if (atmp > ainvnorm) | |
| 1893 ainvnorm = atmp; | |
| 5164 | 1894 } |
| 5681 | 1895 rcond = 1. / ainvnorm / anorm; |
| 5164 | 1896 } |
| 1897 } | |
| 1898 else | |
| 1899 { | |
| 5630 | 1900 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 5164 | 1901 |
| 5275 | 1902 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 1903 { |
| 5630 | 1904 for (octave_idx_type i = 0; i < nm; i++) |
| 5164 | 1905 work[i] = 0.; |
| 5275 | 1906 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) |
| 5164 | 1907 work[b.ridx(i)] = b.data(i); |
| 1908 | |
| 5630 | 1909 for (octave_idx_type k = nc-1; k >= 0; k--) |
| 5164 | 1910 { |
| 1911 if (work[k] != 0.) | |
| 1912 { | |
| 5681 | 1913 if (ridx(cidx(k+1)-1) != k || |
| 1914 data(cidx(k+1)-1) == 0.) | |
| 5164 | 1915 { |
| 1916 err = -2; | |
| 1917 goto triangular_error; | |
| 1918 } | |
| 1919 | |
| 1920 Complex tmp = work[k] / data(cidx(k+1)-1); | |
| 1921 work[k] = tmp; | |
| 5275 | 1922 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) |
| 5164 | 1923 { |
| 5275 | 1924 octave_idx_type iidx = ridx(i); |
| 5164 | 1925 work[iidx] = work[iidx] - tmp * data(i); |
| 1926 } | |
| 1927 } | |
| 1928 } | |
| 1929 | |
| 1930 // Count non-zeros in work vector and adjust space in | |
| 1931 // retval if needed | |
| 5275 | 1932 octave_idx_type new_nnz = 0; |
| 5630 | 1933 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 1934 if (work[i] != 0.) |
| 1935 new_nnz++; | |
| 1936 | |
| 1937 if (ii + new_nnz > x_nz) | |
| 1938 { | |
| 1939 // Resize the sparse matrix | |
| 5275 | 1940 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; |
| 5164 | 1941 retval.change_capacity (sz); |
| 1942 x_nz = sz; | |
| 1943 } | |
| 1944 | |
| 5630 | 1945 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 1946 if (work[i] != 0.) |
| 1947 { | |
| 1948 retval.xridx(ii) = i; | |
| 1949 retval.xdata(ii++) = work[i]; | |
| 1950 } | |
| 1951 retval.xcidx(j+1) = ii; | |
| 1952 } | |
| 1953 | |
| 1954 retval.maybe_compress (); | |
| 1955 | |
| 5681 | 1956 if (calc_cond) |
| 1957 { | |
| 1958 // Calculation of 1-norm of inv(*this) | |
| 1959 for (octave_idx_type i = 0; i < nm; i++) | |
| 1960 work[i] = 0.; | |
| 1961 | |
| 1962 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 1963 { |
| 5681 | 1964 work[j] = 1.; |
| 1965 | |
| 1966 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 1967 { |
| 5681 | 1968 if (work[k] != 0.) |
| 5164 | 1969 { |
| 5681 | 1970 Complex tmp = work[k] / data(cidx(k+1)-1); |
| 1971 work[k] = tmp; | |
| 1972 for (octave_idx_type i = cidx(k); | |
| 1973 i < cidx(k+1)-1; i++) | |
| 1974 { | |
| 1975 octave_idx_type iidx = ridx(i); | |
| 1976 work[iidx] = work[iidx] - tmp * data(i); | |
| 1977 } | |
| 5164 | 1978 } |
| 1979 } | |
| 5681 | 1980 double atmp = 0; |
| 1981 for (octave_idx_type i = 0; i < j+1; i++) | |
| 1982 { | |
| 1983 atmp += std::abs(work[i]); | |
| 1984 work[i] = 0.; | |
| 1985 } | |
| 1986 if (atmp > ainvnorm) | |
| 1987 ainvnorm = atmp; | |
| 5164 | 1988 } |
| 5681 | 1989 rcond = 1. / ainvnorm / anorm; |
| 1990 } | |
| 1991 } | |
| 5164 | 1992 |
| 1993 triangular_error: | |
| 1994 if (err != 0) | |
| 1995 { | |
| 1996 if (sing_handler) | |
| 5681 | 1997 { |
| 1998 sing_handler (rcond); | |
| 1999 mattype.mark_as_rectangular (); | |
| 2000 } | |
| 5164 | 2001 else |
| 2002 (*current_liboctave_error_handler) | |
| 2003 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", | |
| 2004 rcond); | |
| 2005 } | |
| 2006 | |
| 2007 volatile double rcond_plus_one = rcond + 1.0; | |
| 2008 | |
| 2009 if (rcond_plus_one == 1.0 || xisnan (rcond)) | |
| 2010 { | |
| 2011 err = -2; | |
| 2012 | |
| 2013 if (sing_handler) | |
| 5681 | 2014 { |
| 2015 sing_handler (rcond); | |
| 2016 mattype.mark_as_rectangular (); | |
| 2017 } | |
| 5164 | 2018 else |
| 2019 (*current_liboctave_error_handler) | |
| 2020 ("matrix singular to machine precision, rcond = %g", | |
| 2021 rcond); | |
| 2022 } | |
| 2023 } | |
| 2024 else | |
| 2025 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 2026 } | |
| 2027 return retval; | |
| 2028 } | |
| 2029 | |
| 2030 ComplexMatrix | |
| 5785 | 2031 SparseComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b, |
| 5630 | 2032 octave_idx_type& err, double& rcond, |
| 5681 | 2033 solve_singularity_handler sing_handler, |
| 2034 bool calc_cond) const | |
| 5164 | 2035 { |
| 2036 ComplexMatrix retval; | |
| 2037 | |
| 5275 | 2038 octave_idx_type nr = rows (); |
| 2039 octave_idx_type nc = cols (); | |
| 5630 | 2040 octave_idx_type nm = (nc > nr ? nc : nr); |
| 5164 | 2041 err = 0; |
| 2042 | |
| 6924 | 2043 if (nr != b.rows ()) |
| 5164 | 2044 (*current_liboctave_error_handler) |
| 2045 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 2046 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 2047 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
| 5164 | 2048 else |
| 2049 { | |
| 2050 // Print spparms("spumoni") info if requested | |
| 2051 int typ = mattype.type (); | |
| 2052 mattype.info (); | |
| 2053 | |
| 5785 | 2054 if (typ == MatrixType::Permuted_Upper || |
| 2055 typ == MatrixType::Upper) | |
| 5164 | 2056 { |
| 2057 double anorm = 0.; | |
| 2058 double ainvnorm = 0.; | |
| 5275 | 2059 octave_idx_type b_nc = b.cols (); |
| 5681 | 2060 rcond = 1.; |
| 2061 | |
| 2062 if (calc_cond) | |
| 2063 { | |
| 2064 // Calculate the 1-norm of matrix for rcond calculation | |
| 2065 for (octave_idx_type j = 0; j < nc; j++) | |
| 2066 { | |
| 2067 double atmp = 0.; | |
| 2068 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 2069 atmp += std::abs(data(i)); | |
| 2070 if (atmp > anorm) | |
| 2071 anorm = atmp; | |
| 2072 } | |
| 5164 | 2073 } |
| 2074 | |
| 5785 | 2075 if (typ == MatrixType::Permuted_Upper) |
| 5164 | 2076 { |
| 5630 | 2077 retval.resize (nc, b_nc); |
| 5322 | 2078 octave_idx_type *perm = mattype.triangular_perm (); |
| 5630 | 2079 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 5164 | 2080 |
| 5275 | 2081 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 2082 { |
| 5275 | 2083 for (octave_idx_type i = 0; i < nr; i++) |
| 5164 | 2084 work[i] = b(i,j); |
| 5630 | 2085 for (octave_idx_type i = nr; i < nc; i++) |
| 2086 work[i] = 0.; | |
| 2087 | |
| 2088 for (octave_idx_type k = nc-1; k >= 0; k--) | |
| 5164 | 2089 { |
| 5322 | 2090 octave_idx_type kidx = perm[k]; |
| 2091 | |
| 2092 if (work[k] != 0.) | |
| 5164 | 2093 { |
| 5681 | 2094 if (ridx(cidx(kidx+1)-1) != k || |
| 2095 data(cidx(kidx+1)-1) == 0.) | |
| 5164 | 2096 { |
| 2097 err = -2; | |
| 2098 goto triangular_error; | |
| 2099 } | |
| 2100 | |
| 5322 | 2101 Complex tmp = work[k] / data(cidx(kidx+1)-1); |
| 2102 work[k] = tmp; | |
| 2103 for (octave_idx_type i = cidx(kidx); | |
| 2104 i < cidx(kidx+1)-1; i++) | |
| 5164 | 2105 { |
| 5322 | 2106 octave_idx_type iidx = ridx(i); |
| 2107 work[iidx] = work[iidx] - tmp * data(i); | |
| 5164 | 2108 } |
| 2109 } | |
| 2110 } | |
| 2111 | |
| 5630 | 2112 for (octave_idx_type i = 0; i < nc; i++) |
| 5322 | 2113 retval (perm[i], j) = work[i]; |
| 5164 | 2114 } |
| 2115 | |
| 5681 | 2116 if (calc_cond) |
| 2117 { | |
| 2118 // Calculation of 1-norm of inv(*this) | |
| 2119 for (octave_idx_type i = 0; i < nm; i++) | |
| 2120 work[i] = 0.; | |
| 2121 | |
| 2122 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 2123 { |
| 5681 | 2124 work[j] = 1.; |
| 2125 | |
| 2126 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 2127 { |
| 5681 | 2128 octave_idx_type iidx = perm[k]; |
| 2129 | |
| 2130 if (work[k] != 0.) | |
| 5164 | 2131 { |
| 5681 | 2132 Complex tmp = work[k] / data(cidx(iidx+1)-1); |
| 2133 work[k] = tmp; | |
| 2134 for (octave_idx_type i = cidx(iidx); | |
| 2135 i < cidx(iidx+1)-1; i++) | |
| 2136 { | |
| 2137 octave_idx_type idx2 = ridx(i); | |
| 2138 work[idx2] = work[idx2] - tmp * data(i); | |
| 2139 } | |
| 5164 | 2140 } |
| 2141 } | |
| 5681 | 2142 double atmp = 0; |
| 2143 for (octave_idx_type i = 0; i < j+1; i++) | |
| 2144 { | |
| 2145 atmp += std::abs(work[i]); | |
| 2146 work[i] = 0.; | |
| 2147 } | |
| 2148 if (atmp > ainvnorm) | |
| 2149 ainvnorm = atmp; | |
| 5164 | 2150 } |
| 5681 | 2151 rcond = 1. / ainvnorm / anorm; |
| 5164 | 2152 } |
| 2153 } | |
| 2154 else | |
| 2155 { | |
| 5630 | 2156 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 2157 retval.resize (nc, b_nc); | |
| 5164 | 2158 |
| 5275 | 2159 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 2160 { |
| 5630 | 2161 for (octave_idx_type i = 0; i < nr; i++) |
| 2162 work[i] = b(i,j); | |
| 2163 for (octave_idx_type i = nr; i < nc; i++) | |
| 2164 work[i] = 0.; | |
| 2165 | |
| 2166 for (octave_idx_type k = nc-1; k >= 0; k--) | |
| 5164 | 2167 { |
| 5630 | 2168 if (work[k] != 0.) |
| 5164 | 2169 { |
| 5681 | 2170 if (ridx(cidx(k+1)-1) != k || |
| 2171 data(cidx(k+1)-1) == 0.) | |
| 5164 | 2172 { |
| 2173 err = -2; | |
| 2174 goto triangular_error; | |
| 2175 } | |
| 2176 | |
| 5630 | 2177 Complex tmp = work[k] / data(cidx(k+1)-1); |
| 2178 work[k] = tmp; | |
| 5275 | 2179 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) |
| 5164 | 2180 { |
| 5275 | 2181 octave_idx_type iidx = ridx(i); |
| 5630 | 2182 work[iidx] = work[iidx] - tmp * data(i); |
| 5164 | 2183 } |
| 2184 } | |
| 2185 } | |
| 5630 | 2186 |
| 2187 for (octave_idx_type i = 0; i < nc; i++) | |
| 2188 retval.xelem (i, j) = work[i]; | |
| 5164 | 2189 } |
| 2190 | |
| 5681 | 2191 if (calc_cond) |
| 2192 { | |
| 2193 // Calculation of 1-norm of inv(*this) | |
| 2194 for (octave_idx_type i = 0; i < nm; i++) | |
| 2195 work[i] = 0.; | |
| 2196 | |
| 2197 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 2198 { |
| 5681 | 2199 work[j] = 1.; |
| 2200 | |
| 2201 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 2202 { |
| 5681 | 2203 if (work[k] != 0.) |
| 5164 | 2204 { |
| 5681 | 2205 Complex tmp = work[k] / data(cidx(k+1)-1); |
| 2206 work[k] = tmp; | |
| 2207 for (octave_idx_type i = cidx(k); | |
| 2208 i < cidx(k+1)-1; i++) | |
| 2209 { | |
| 2210 octave_idx_type iidx = ridx(i); | |
| 2211 work[iidx] = work[iidx] - tmp * data(i); | |
| 2212 } | |
| 5164 | 2213 } |
| 2214 } | |
| 5681 | 2215 double atmp = 0; |
| 2216 for (octave_idx_type i = 0; i < j+1; i++) | |
| 2217 { | |
| 2218 atmp += std::abs(work[i]); | |
| 2219 work[i] = 0.; | |
| 2220 } | |
| 2221 if (atmp > ainvnorm) | |
| 2222 ainvnorm = atmp; | |
| 5164 | 2223 } |
| 5681 | 2224 rcond = 1. / ainvnorm / anorm; |
| 2225 } | |
| 2226 } | |
| 5164 | 2227 |
| 2228 triangular_error: | |
| 2229 if (err != 0) | |
| 2230 { | |
| 2231 if (sing_handler) | |
| 5681 | 2232 { |
| 2233 sing_handler (rcond); | |
| 2234 mattype.mark_as_rectangular (); | |
| 2235 } | |
| 5164 | 2236 else |
| 2237 (*current_liboctave_error_handler) | |
| 2238 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", | |
| 2239 rcond); | |
| 2240 } | |
| 2241 | |
| 2242 volatile double rcond_plus_one = rcond + 1.0; | |
| 2243 | |
| 2244 if (rcond_plus_one == 1.0 || xisnan (rcond)) | |
| 2245 { | |
| 2246 err = -2; | |
| 2247 | |
| 2248 if (sing_handler) | |
| 5681 | 2249 { |
| 2250 sing_handler (rcond); | |
| 2251 mattype.mark_as_rectangular (); | |
| 2252 } | |
| 5164 | 2253 else |
| 2254 (*current_liboctave_error_handler) | |
| 2255 ("matrix singular to machine precision, rcond = %g", | |
| 2256 rcond); | |
| 2257 } | |
| 2258 } | |
| 2259 else | |
| 2260 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 2261 } | |
| 2262 | |
| 2263 return retval; | |
| 2264 } | |
| 2265 | |
| 2266 SparseComplexMatrix | |
| 5785 | 2267 SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseComplexMatrix& b, |
| 5630 | 2268 octave_idx_type& err, double& rcond, |
| 5681 | 2269 solve_singularity_handler sing_handler, |
| 2270 bool calc_cond) const | |
| 5164 | 2271 { |
| 2272 SparseComplexMatrix retval; | |
| 2273 | |
| 5275 | 2274 octave_idx_type nr = rows (); |
| 2275 octave_idx_type nc = cols (); | |
| 5630 | 2276 octave_idx_type nm = (nc > nr ? nc : nr); |
| 5164 | 2277 err = 0; |
| 2278 | |
| 6924 | 2279 if (nr != b.rows ()) |
| 5164 | 2280 (*current_liboctave_error_handler) |
| 2281 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 2282 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 2283 retval = SparseComplexMatrix (nc, b.cols ()); | |
| 5164 | 2284 else |
| 2285 { | |
| 2286 // Print spparms("spumoni") info if requested | |
| 2287 int typ = mattype.type (); | |
| 2288 mattype.info (); | |
| 2289 | |
| 5785 | 2290 if (typ == MatrixType::Permuted_Upper || |
| 2291 typ == MatrixType::Upper) | |
| 5164 | 2292 { |
| 2293 double anorm = 0.; | |
| 2294 double ainvnorm = 0.; | |
| 5681 | 2295 rcond = 1.; |
| 2296 | |
| 2297 if (calc_cond) | |
| 2298 { | |
| 2299 // Calculate the 1-norm of matrix for rcond calculation | |
| 2300 for (octave_idx_type j = 0; j < nc; j++) | |
| 2301 { | |
| 2302 double atmp = 0.; | |
| 2303 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 2304 atmp += std::abs(data(i)); | |
| 2305 if (atmp > anorm) | |
| 2306 anorm = atmp; | |
| 2307 } | |
| 5164 | 2308 } |
| 2309 | |
| 5275 | 2310 octave_idx_type b_nc = b.cols (); |
| 5681 | 2311 octave_idx_type b_nz = b.nnz (); |
| 5630 | 2312 retval = SparseComplexMatrix (nc, b_nc, b_nz); |
| 5164 | 2313 retval.xcidx(0) = 0; |
| 5275 | 2314 octave_idx_type ii = 0; |
| 2315 octave_idx_type x_nz = b_nz; | |
| 5164 | 2316 |
| 5785 | 2317 if (typ == MatrixType::Permuted_Upper) |
| 5164 | 2318 { |
| 5322 | 2319 octave_idx_type *perm = mattype.triangular_perm (); |
| 5630 | 2320 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 2321 | |
| 2322 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); | |
| 2323 for (octave_idx_type i = 0; i < nc; i++) | |
| 5322 | 2324 rperm[perm[i]] = i; |
| 5164 | 2325 |
| 5275 | 2326 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 2327 { |
| 5630 | 2328 for (octave_idx_type i = 0; i < nm; i++) |
| 5164 | 2329 work[i] = 0.; |
| 5275 | 2330 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) |
| 5164 | 2331 work[b.ridx(i)] = b.data(i); |
| 2332 | |
| 5630 | 2333 for (octave_idx_type k = nc-1; k >= 0; k--) |
| 5164 | 2334 { |
| 5322 | 2335 octave_idx_type kidx = perm[k]; |
| 2336 | |
| 2337 if (work[k] != 0.) | |
| 5164 | 2338 { |
| 5681 | 2339 if (ridx(cidx(kidx+1)-1) != k || |
| 2340 data(cidx(kidx+1)-1) == 0.) | |
| 5164 | 2341 { |
| 2342 err = -2; | |
| 2343 goto triangular_error; | |
| 2344 } | |
| 2345 | |
| 5322 | 2346 Complex tmp = work[k] / data(cidx(kidx+1)-1); |
| 2347 work[k] = tmp; | |
| 2348 for (octave_idx_type i = cidx(kidx); | |
| 2349 i < cidx(kidx+1)-1; i++) | |
| 5164 | 2350 { |
| 5322 | 2351 octave_idx_type iidx = ridx(i); |
| 2352 work[iidx] = work[iidx] - tmp * data(i); | |
| 5164 | 2353 } |
| 2354 } | |
| 2355 } | |
| 2356 | |
| 2357 // Count non-zeros in work vector and adjust space in | |
| 2358 // retval if needed | |
| 5275 | 2359 octave_idx_type new_nnz = 0; |
| 5630 | 2360 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 2361 if (work[i] != 0.) |
| 2362 new_nnz++; | |
| 2363 | |
| 2364 if (ii + new_nnz > x_nz) | |
| 2365 { | |
| 2366 // Resize the sparse matrix | |
| 5275 | 2367 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; |
| 5164 | 2368 retval.change_capacity (sz); |
| 2369 x_nz = sz; | |
| 2370 } | |
| 2371 | |
| 5630 | 2372 for (octave_idx_type i = 0; i < nc; i++) |
| 5322 | 2373 if (work[rperm[i]] != 0.) |
| 5164 | 2374 { |
| 2375 retval.xridx(ii) = i; | |
| 5322 | 2376 retval.xdata(ii++) = work[rperm[i]]; |
| 5164 | 2377 } |
| 2378 retval.xcidx(j+1) = ii; | |
| 2379 } | |
| 2380 | |
| 2381 retval.maybe_compress (); | |
| 2382 | |
| 5681 | 2383 if (calc_cond) |
| 2384 { | |
| 2385 // Calculation of 1-norm of inv(*this) | |
| 2386 for (octave_idx_type i = 0; i < nm; i++) | |
| 2387 work[i] = 0.; | |
| 2388 | |
| 2389 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 2390 { |
| 5681 | 2391 work[j] = 1.; |
| 2392 | |
| 2393 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 2394 { |
| 5681 | 2395 octave_idx_type iidx = perm[k]; |
| 2396 | |
| 2397 if (work[k] != 0.) | |
| 5164 | 2398 { |
| 5681 | 2399 Complex tmp = work[k] / data(cidx(iidx+1)-1); |
| 2400 work[k] = tmp; | |
| 2401 for (octave_idx_type i = cidx(iidx); | |
| 2402 i < cidx(iidx+1)-1; i++) | |
| 2403 { | |
| 2404 octave_idx_type idx2 = ridx(i); | |
| 2405 work[idx2] = work[idx2] - tmp * data(i); | |
| 2406 } | |
| 5164 | 2407 } |
| 2408 } | |
| 5681 | 2409 double atmp = 0; |
| 2410 for (octave_idx_type i = 0; i < j+1; i++) | |
| 2411 { | |
| 2412 atmp += std::abs(work[i]); | |
| 2413 work[i] = 0.; | |
| 2414 } | |
| 2415 if (atmp > ainvnorm) | |
| 2416 ainvnorm = atmp; | |
| 5164 | 2417 } |
| 5681 | 2418 rcond = 1. / ainvnorm / anorm; |
| 5164 | 2419 } |
| 2420 } | |
| 2421 else | |
| 2422 { | |
| 5630 | 2423 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 5164 | 2424 |
| 5275 | 2425 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 2426 { |
| 5630 | 2427 for (octave_idx_type i = 0; i < nm; i++) |
| 5164 | 2428 work[i] = 0.; |
| 5275 | 2429 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) |
| 5164 | 2430 work[b.ridx(i)] = b.data(i); |
| 2431 | |
| 5275 | 2432 for (octave_idx_type k = nr-1; k >= 0; k--) |
| 5164 | 2433 { |
| 2434 if (work[k] != 0.) | |
| 2435 { | |
| 5681 | 2436 if (ridx(cidx(k+1)-1) != k || |
| 2437 data(cidx(k+1)-1) == 0.) | |
| 5164 | 2438 { |
| 2439 err = -2; | |
| 2440 goto triangular_error; | |
| 2441 } | |
| 2442 | |
| 2443 Complex tmp = work[k] / data(cidx(k+1)-1); | |
| 2444 work[k] = tmp; | |
| 5275 | 2445 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) |
| 5164 | 2446 { |
| 5275 | 2447 octave_idx_type iidx = ridx(i); |
| 5164 | 2448 work[iidx] = work[iidx] - tmp * data(i); |
| 2449 } | |
| 2450 } | |
| 2451 } | |
| 2452 | |
| 2453 // Count non-zeros in work vector and adjust space in | |
| 2454 // retval if needed | |
| 5275 | 2455 octave_idx_type new_nnz = 0; |
| 5630 | 2456 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 2457 if (work[i] != 0.) |
| 2458 new_nnz++; | |
| 2459 | |
| 2460 if (ii + new_nnz > x_nz) | |
| 2461 { | |
| 2462 // Resize the sparse matrix | |
| 5275 | 2463 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; |
| 5164 | 2464 retval.change_capacity (sz); |
| 2465 x_nz = sz; | |
| 2466 } | |
| 2467 | |
| 5630 | 2468 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 2469 if (work[i] != 0.) |
| 2470 { | |
| 2471 retval.xridx(ii) = i; | |
| 2472 retval.xdata(ii++) = work[i]; | |
| 2473 } | |
| 2474 retval.xcidx(j+1) = ii; | |
| 2475 } | |
| 2476 | |
| 2477 retval.maybe_compress (); | |
| 2478 | |
| 5681 | 2479 if (calc_cond) |
| 2480 { | |
| 2481 // Calculation of 1-norm of inv(*this) | |
| 2482 for (octave_idx_type i = 0; i < nm; i++) | |
| 2483 work[i] = 0.; | |
| 2484 | |
| 2485 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 2486 { |
| 5681 | 2487 work[j] = 1.; |
| 2488 | |
| 2489 for (octave_idx_type k = j; k >= 0; k--) | |
| 5164 | 2490 { |
| 5681 | 2491 if (work[k] != 0.) |
| 5164 | 2492 { |
| 5681 | 2493 Complex tmp = work[k] / data(cidx(k+1)-1); |
| 2494 work[k] = tmp; | |
| 2495 for (octave_idx_type i = cidx(k); | |
| 2496 i < cidx(k+1)-1; i++) | |
| 2497 { | |
| 2498 octave_idx_type iidx = ridx(i); | |
| 2499 work[iidx] = work[iidx] - tmp * data(i); | |
| 2500 } | |
| 5164 | 2501 } |
| 2502 } | |
| 5681 | 2503 double atmp = 0; |
| 2504 for (octave_idx_type i = 0; i < j+1; i++) | |
| 2505 { | |
| 2506 atmp += std::abs(work[i]); | |
| 2507 work[i] = 0.; | |
| 2508 } | |
| 2509 if (atmp > ainvnorm) | |
| 2510 ainvnorm = atmp; | |
| 5164 | 2511 } |
| 5681 | 2512 rcond = 1. / ainvnorm / anorm; |
| 2513 } | |
| 2514 } | |
| 5164 | 2515 |
| 2516 triangular_error: | |
| 2517 if (err != 0) | |
| 2518 { | |
| 2519 if (sing_handler) | |
| 5681 | 2520 { |
| 2521 sing_handler (rcond); | |
| 2522 mattype.mark_as_rectangular (); | |
| 2523 } | |
| 5164 | 2524 else |
| 2525 (*current_liboctave_error_handler) | |
| 2526 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", | |
| 2527 rcond); | |
| 2528 } | |
| 2529 | |
| 2530 volatile double rcond_plus_one = rcond + 1.0; | |
| 2531 | |
| 2532 if (rcond_plus_one == 1.0 || xisnan (rcond)) | |
| 2533 { | |
| 2534 err = -2; | |
| 2535 | |
| 2536 if (sing_handler) | |
| 5681 | 2537 { |
| 2538 sing_handler (rcond); | |
| 2539 mattype.mark_as_rectangular (); | |
| 2540 } | |
| 5164 | 2541 else |
| 2542 (*current_liboctave_error_handler) | |
| 2543 ("matrix singular to machine precision, rcond = %g", | |
| 2544 rcond); | |
| 2545 } | |
| 2546 } | |
| 2547 else | |
| 2548 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 2549 } | |
| 2550 | |
| 2551 return retval; | |
| 2552 } | |
| 2553 | |
| 2554 ComplexMatrix | |
| 5785 | 2555 SparseComplexMatrix::ltsolve (MatrixType &mattype, const Matrix& b, |
| 5630 | 2556 octave_idx_type& err, double& rcond, |
| 5681 | 2557 solve_singularity_handler sing_handler, |
| 2558 bool calc_cond) const | |
| 5164 | 2559 { |
| 2560 ComplexMatrix retval; | |
| 2561 | |
| 5275 | 2562 octave_idx_type nr = rows (); |
| 2563 octave_idx_type nc = cols (); | |
| 5630 | 2564 octave_idx_type nm = (nc > nr ? nc : nr); |
| 5164 | 2565 err = 0; |
| 2566 | |
| 6924 | 2567 if (nr != b.rows ()) |
| 5164 | 2568 (*current_liboctave_error_handler) |
| 2569 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 2570 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 2571 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
| 5164 | 2572 else |
| 2573 { | |
| 2574 // Print spparms("spumoni") info if requested | |
| 2575 int typ = mattype.type (); | |
| 2576 mattype.info (); | |
| 2577 | |
| 5785 | 2578 if (typ == MatrixType::Permuted_Lower || |
| 2579 typ == MatrixType::Lower) | |
| 5164 | 2580 { |
| 2581 double anorm = 0.; | |
| 2582 double ainvnorm = 0.; | |
| 5630 | 2583 octave_idx_type b_nc = b.cols (); |
| 5681 | 2584 rcond = 1.; |
| 2585 | |
| 2586 if (calc_cond) | |
| 2587 { | |
| 2588 // Calculate the 1-norm of matrix for rcond calculation | |
| 2589 for (octave_idx_type j = 0; j < nc; j++) | |
| 2590 { | |
| 2591 double atmp = 0.; | |
| 2592 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 2593 atmp += std::abs(data(i)); | |
| 2594 if (atmp > anorm) | |
| 2595 anorm = atmp; | |
| 2596 } | |
| 5164 | 2597 } |
| 2598 | |
| 5785 | 2599 if (typ == MatrixType::Permuted_Lower) |
| 5164 | 2600 { |
| 5630 | 2601 retval.resize (nc, b_nc); |
| 2602 OCTAVE_LOCAL_BUFFER (Complex, work, nm); | |
| 5322 | 2603 octave_idx_type *perm = mattype.triangular_perm (); |
| 5164 | 2604 |
| 5630 | 2605 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 2606 { |
| 5630 | 2607 for (octave_idx_type i = 0; i < nm; i++) |
| 2608 work[i] = 0.; | |
| 5275 | 2609 for (octave_idx_type i = 0; i < nr; i++) |
| 5322 | 2610 work[perm[i]] = b(i,j); |
| 5164 | 2611 |
| 5630 | 2612 for (octave_idx_type k = 0; k < nc; k++) |
| 5164 | 2613 { |
| 5322 | 2614 if (work[k] != 0.) |
| 5164 | 2615 { |
| 5322 | 2616 octave_idx_type minr = nr; |
| 2617 octave_idx_type mini = 0; | |
| 2618 | |
| 2619 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) | |
| 2620 if (perm[ridx(i)] < minr) | |
| 2621 { | |
| 2622 minr = perm[ridx(i)]; | |
| 2623 mini = i; | |
| 2624 } | |
| 2625 | |
| 5681 | 2626 if (minr != k || data (mini) == 0.) |
| 5164 | 2627 { |
| 2628 err = -2; | |
| 2629 goto triangular_error; | |
| 2630 } | |
| 2631 | |
| 5322 | 2632 Complex tmp = work[k] / data(mini); |
| 2633 work[k] = tmp; | |
| 2634 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) | |
| 5164 | 2635 { |
| 5322 | 2636 if (i == mini) |
| 2637 continue; | |
| 2638 | |
| 2639 octave_idx_type iidx = perm[ridx(i)]; | |
| 2640 work[iidx] = work[iidx] - tmp * data(i); | |
| 5164 | 2641 } |
| 2642 } | |
| 2643 } | |
| 2644 | |
| 5630 | 2645 for (octave_idx_type i = 0; i < nc; i++) |
| 5322 | 2646 retval (i, j) = work[i]; |
| 5164 | 2647 } |
| 2648 | |
| 5681 | 2649 if (calc_cond) |
| 2650 { | |
| 2651 // Calculation of 1-norm of inv(*this) | |
| 2652 for (octave_idx_type i = 0; i < nm; i++) | |
| 2653 work[i] = 0.; | |
| 2654 | |
| 2655 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 2656 { |
| 5681 | 2657 work[j] = 1.; |
| 2658 | |
| 2659 for (octave_idx_type k = 0; k < nc; k++) | |
| 5164 | 2660 { |
| 5681 | 2661 if (work[k] != 0.) |
| 5164 | 2662 { |
| 5681 | 2663 octave_idx_type minr = nr; |
| 2664 octave_idx_type mini = 0; | |
| 2665 | |
| 2666 for (octave_idx_type i = cidx(k); | |
| 2667 i < cidx(k+1); i++) | |
| 2668 if (perm[ridx(i)] < minr) | |
| 2669 { | |
| 2670 minr = perm[ridx(i)]; | |
| 2671 mini = i; | |
| 2672 } | |
| 2673 | |
| 2674 Complex tmp = work[k] / data(mini); | |
| 2675 work[k] = tmp; | |
| 2676 for (octave_idx_type i = cidx(k); | |
| 2677 i < cidx(k+1); i++) | |
| 2678 { | |
| 2679 if (i == mini) | |
| 2680 continue; | |
| 2681 | |
| 2682 octave_idx_type iidx = perm[ridx(i)]; | |
| 2683 work[iidx] = work[iidx] - tmp * data(i); | |
| 2684 } | |
| 5164 | 2685 } |
| 2686 } | |
| 5681 | 2687 |
| 2688 double atmp = 0; | |
| 2689 for (octave_idx_type i = j; i < nc; i++) | |
| 2690 { | |
| 2691 atmp += std::abs(work[i]); | |
| 2692 work[i] = 0.; | |
| 2693 } | |
| 2694 if (atmp > ainvnorm) | |
| 2695 ainvnorm = atmp; | |
| 5164 | 2696 } |
| 5681 | 2697 rcond = 1. / ainvnorm / anorm; |
| 5164 | 2698 } |
| 2699 } | |
| 2700 else | |
| 2701 { | |
| 5630 | 2702 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 2703 retval.resize (nc, b_nc, 0.); | |
| 2704 | |
| 2705 for (octave_idx_type j = 0; j < b_nc; j++) | |
| 5164 | 2706 { |
| 5630 | 2707 for (octave_idx_type i = 0; i < nr; i++) |
| 2708 work[i] = b(i,j); | |
| 2709 for (octave_idx_type i = nr; i < nc; i++) | |
| 2710 work[i] = 0.; | |
| 2711 for (octave_idx_type k = 0; k < nc; k++) | |
| 5164 | 2712 { |
| 5630 | 2713 if (work[k] != 0.) |
| 5164 | 2714 { |
| 5681 | 2715 if (ridx(cidx(k)) != k || |
| 2716 data(cidx(k)) == 0.) | |
| 5164 | 2717 { |
| 2718 err = -2; | |
| 2719 goto triangular_error; | |
| 2720 } | |
| 2721 | |
| 5630 | 2722 Complex tmp = work[k] / data(cidx(k)); |
| 2723 work[k] = tmp; | |
| 5275 | 2724 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) |
| 5164 | 2725 { |
| 5275 | 2726 octave_idx_type iidx = ridx(i); |
| 5630 | 2727 work[iidx] = work[iidx] - tmp * data(i); |
| 5164 | 2728 } |
| 2729 } | |
| 2730 } | |
| 5630 | 2731 for (octave_idx_type i = 0; i < nc; i++) |
| 2732 retval.xelem (i, j) = work[i]; | |
| 5164 | 2733 } |
| 2734 | |
| 5681 | 2735 if (calc_cond) |
| 2736 { | |
| 2737 // Calculation of 1-norm of inv(*this) | |
| 2738 for (octave_idx_type i = 0; i < nm; i++) | |
| 2739 work[i] = 0.; | |
| 2740 | |
| 2741 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 2742 { |
| 5681 | 2743 work[j] = 1.; |
| 2744 | |
| 2745 for (octave_idx_type k = j; k < nc; k++) | |
| 5164 | 2746 { |
| 5681 | 2747 |
| 2748 if (work[k] != 0.) | |
| 5164 | 2749 { |
| 5681 | 2750 Complex tmp = work[k] / data(cidx(k)); |
| 2751 work[k] = tmp; | |
| 2752 for (octave_idx_type i = cidx(k)+1; | |
| 2753 i < cidx(k+1); i++) | |
| 2754 { | |
| 2755 octave_idx_type iidx = ridx(i); | |
| 2756 work[iidx] = work[iidx] - tmp * data(i); | |
| 2757 } | |
| 5164 | 2758 } |
| 2759 } | |
| 5681 | 2760 double atmp = 0; |
| 2761 for (octave_idx_type i = j; i < nc; i++) | |
| 2762 { | |
| 2763 atmp += std::abs(work[i]); | |
| 2764 work[i] = 0.; | |
| 2765 } | |
| 2766 if (atmp > ainvnorm) | |
| 2767 ainvnorm = atmp; | |
| 5164 | 2768 } |
| 5681 | 2769 rcond = 1. / ainvnorm / anorm; |
| 2770 } | |
| 2771 } | |
| 5164 | 2772 triangular_error: |
| 2773 if (err != 0) | |
| 2774 { | |
| 2775 if (sing_handler) | |
| 5681 | 2776 { |
| 2777 sing_handler (rcond); | |
| 2778 mattype.mark_as_rectangular (); | |
| 2779 } | |
| 5164 | 2780 else |
| 2781 (*current_liboctave_error_handler) | |
| 2782 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", | |
| 2783 rcond); | |
| 2784 } | |
| 2785 | |
| 2786 volatile double rcond_plus_one = rcond + 1.0; | |
| 2787 | |
| 2788 if (rcond_plus_one == 1.0 || xisnan (rcond)) | |
| 2789 { | |
| 2790 err = -2; | |
| 2791 | |
| 2792 if (sing_handler) | |
| 5681 | 2793 { |
| 2794 sing_handler (rcond); | |
| 2795 mattype.mark_as_rectangular (); | |
| 2796 } | |
| 5164 | 2797 else |
| 2798 (*current_liboctave_error_handler) | |
| 2799 ("matrix singular to machine precision, rcond = %g", | |
| 2800 rcond); | |
| 2801 } | |
| 2802 } | |
| 2803 else | |
| 2804 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
| 2805 } | |
| 2806 | |
| 2807 return retval; | |
| 2808 } | |
| 2809 | |
| 2810 SparseComplexMatrix | |
| 5785 | 2811 SparseComplexMatrix::ltsolve (MatrixType &mattype, const SparseMatrix& b, |
| 5630 | 2812 octave_idx_type& err, double& rcond, |
| 5681 | 2813 solve_singularity_handler sing_handler, |
| 2814 bool calc_cond) const | |
| 5164 | 2815 { |
| 2816 SparseComplexMatrix retval; | |
| 2817 | |
| 5275 | 2818 octave_idx_type nr = rows (); |
| 2819 octave_idx_type nc = cols (); | |
| 5630 | 2820 octave_idx_type nm = (nc > nr ? nc : nr); |
| 2821 | |
| 5164 | 2822 err = 0; |
| 2823 | |
| 6924 | 2824 if (nr != b.rows ()) |
| 5164 | 2825 (*current_liboctave_error_handler) |
| 2826 ("matrix dimension mismatch solution of linear equations"); | |
| 6924 | 2827 else if (nr == 0 || nc == 0 || b.cols () == 0) |
| 2828 retval = SparseComplexMatrix (nc, b.cols ()); | |
| 5164 | 2829 else |
| 2830 { | |
| 2831 // Print spparms("spumoni") info if requested | |
| 2832 int typ = mattype.type (); | |
| 2833 mattype.info (); | |
| 2834 | |
| 5785 | 2835 if (typ == MatrixType::Permuted_Lower || |
| 2836 typ == MatrixType::Lower) | |
| 5164 | 2837 { |
| 2838 double anorm = 0.; | |
| 2839 double ainvnorm = 0.; | |
| 5681 | 2840 rcond = 1.; |
| 2841 | |
| 2842 if (calc_cond) | |
| 2843 { | |
| 2844 // Calculate the 1-norm of matrix for rcond calculation | |
| 2845 for (octave_idx_type j = 0; j < nc; j++) | |
| 2846 { | |
| 2847 double atmp = 0.; | |
| 2848 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) | |
| 2849 atmp += std::abs(data(i)); | |
| 2850 if (atmp > anorm) | |
| 2851 anorm = atmp; | |
| 2852 } | |
| 5164 | 2853 } |
| 2854 | |
| 5275 | 2855 octave_idx_type b_nc = b.cols (); |
| 5681 | 2856 octave_idx_type b_nz = b.nnz (); |
| 5630 | 2857 retval = SparseComplexMatrix (nc, b_nc, b_nz); |
| 5164 | 2858 retval.xcidx(0) = 0; |
| 5275 | 2859 octave_idx_type ii = 0; |
| 2860 octave_idx_type x_nz = b_nz; | |
| 5164 | 2861 |
| 5785 | 2862 if (typ == MatrixType::Permuted_Lower) |
| 5164 | 2863 { |
| 5630 | 2864 OCTAVE_LOCAL_BUFFER (Complex, work, nm); |
| 5322 | 2865 octave_idx_type *perm = mattype.triangular_perm (); |
| 5164 | 2866 |
| 5275 | 2867 for (octave_idx_type j = 0; j < b_nc; j++) |
| 5164 | 2868 { |
| 5630 | 2869 for (octave_idx_type i = 0; i < nm; i++) |
| 5164 | 2870 work[i] = 0.; |
| 5275 | 2871 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) |
| 5322 | 2872 work[perm[b.ridx(i)]] = b.data(i); |
| 5164 | 2873 |
| 5630 | 2874 for (octave_idx_type k = 0; k < nc; k++) |
| 5164 | 2875 { |
| 5322 | 2876 if (work[k] != 0.) |
| 5164 | 2877 { |
| 5322 | 2878 octave_idx_type minr = nr; |
| 2879 octave_idx_type mini = 0; | |
| 2880 | |
| 2881 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) | |
| 2882 if (perm[ridx(i)] < minr) | |
| 2883 { | |
| 2884 minr = perm[ridx(i)]; | |
| 2885 mini = i; | |
| 2886 } | |
| 2887 | |
| 5681 | 2888 if (minr != k || data (mini) == 0.) |
| 5164 | 2889 { |
| 2890 err = -2; | |
| 2891 goto triangular_error; | |
| 2892 } | |
| 2893 | |
| 5322 | 2894 Complex tmp = work[k] / data(mini); |
| 2895 work[k] = tmp; | |
| 2896 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) | |
| 5164 | 2897 { |
| 5322 | 2898 if (i == mini) |
| 2899 continue; | |
| 2900 | |
| 2901 octave_idx_type iidx = perm[ridx(i)]; | |
| 2902 work[iidx] = work[iidx] - tmp * data(i); | |
| 5164 | 2903 } |
| 2904 } | |
| 2905 } | |
| 2906 | |
| 2907 // Count non-zeros in work vector and adjust space in | |
| 2908 // retval if needed | |
| 5275 | 2909 octave_idx_type new_nnz = 0; |
| 5630 | 2910 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 2911 if (work[i] != 0.) |
| 2912 new_nnz++; | |
| 2913 | |
| 2914 if (ii + new_nnz > x_nz) | |
| 2915 { | |
| 2916 // Resize the sparse matrix | |
| 5275 | 2917 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; |
| 5164 | 2918 retval.change_capacity (sz); |
| 2919 x_nz = sz; | |
| 2920 } | |
| 2921 | |
| 5630 | 2922 for (octave_idx_type i = 0; i < nc; i++) |
| 5322 | 2923 if (work[i] != 0.) |
| 5164 | 2924 { |
| 2925 retval.xridx(ii) = i; | |
| 5322 | 2926 retval.xdata(ii++) = work[i]; |
| 5164 | 2927 } |
| 2928 retval.xcidx(j+1) = ii; | |
| 2929 } | |
| 2930 | |
| 2931 retval.maybe_compress (); | |
| 2932 | |
| 5681 | 2933 if (calc_cond) |
| 2934 { | |
| 2935 // Calculation of 1-norm of inv(*this) | |
| 2936 for (octave_idx_type i = 0; i < nm; i++) | |
| 2937 work[i] = 0.; | |
| 2938 | |
| 2939 for (octave_idx_type j = 0; j < nr; j++) | |
| 5164 | 2940 { |
| 5681 | 2941 work[j] = 1.; |
| 2942 | |
| 2943 for (octave_idx_type k = 0; k < nc; k++) | |
| 5164 | 2944 { |
| 5681 | 2945 if (work[k] != 0.) |
| 5164 | 2946 { |
| 5681 | 2947 octave_idx_type minr = nr; |
| 2948 octave_idx_type mini = 0; | |
| 2949 | |
| 2950 for (octave_idx_type i = cidx(k); | |
| 2951 i < cidx(k+1); i++) | |
| 2952 if (perm[ridx(i)] < minr) | |
| 2953 { | |
| 2954 minr = perm[ridx(i)]; | |
| 2955 mini = i; | |
| 2956 } | |
| 2957 | |
| 2958 Complex tmp = work[k] / data(mini); | |
| 2959 work[k] = tmp; | |
| 2960 for (octave_idx_type i = cidx(k); | |
| 2961 i < cidx(k+1); i++) | |
| 2962 { | |
| 2963 if (i == mini) | |
| 2964 continue; | |
| 2965 | |
| 2966 octave_idx_type iidx = perm[ridx(i)]; | |
| 2967 work[iidx] = work[iidx] - tmp * data(i); | |
| 2968 } | |
| 5164 | 2969 } |
| 2970 } | |
| 5681 | 2971 |
| 2972 double atmp = 0; | |
| 2973 for (octave_idx_type i = j; i < nc; i++) | |
| 2974 { | |
| 2975 atmp += std::abs(work[i]); | |
| 2976 work[i] = 0.; | |
| 2977 } | |
|
233d98d95659
[project @ 2006-03-16 17:48:55 by dbateman]
dbateman |