00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifdef HAVE_CONFIG_H
00026 #include <config.h>
00027 #endif
00028
00029 #include <cfloat>
00030
00031 #include <iostream>
00032 #include <vector>
00033
00034 #include "quit.h"
00035 #include "lo-ieee.h"
00036 #include "lo-mappers.h"
00037 #include "f77-fcn.h"
00038 #include "dRowVector.h"
00039 #include "oct-locbuf.h"
00040
00041 #include "dDiagMatrix.h"
00042 #include "CDiagMatrix.h"
00043 #include "CSparse.h"
00044 #include "boolSparse.h"
00045 #include "dSparse.h"
00046 #include "functor.h"
00047 #include "oct-spparms.h"
00048 #include "SparseCmplxLU.h"
00049 #include "oct-sparse.h"
00050 #include "sparse-util.h"
00051 #include "SparseCmplxCHOL.h"
00052 #include "SparseCmplxQR.h"
00053
00054 #include "Sparse-diag-op-defs.h"
00055
00056 #include "Sparse-perm-op-defs.h"
00057 #include "mx-inlines.cc"
00058
00059
00060
00061
00062 #ifndef USE_QRSOLVE
00063 #include "sparse-dmsolve.cc"
00064 #endif
00065
00066
00067 extern "C"
00068 {
00069 F77_RET_T
00070 F77_FUNC (zgbtrf, ZGBTRF) (const octave_idx_type&, const octave_idx_type&,
00071 const octave_idx_type&, const octave_idx_type&,
00072 Complex*, const octave_idx_type&,
00073 octave_idx_type*, octave_idx_type&);
00074
00075 F77_RET_T
00076 F77_FUNC (zgbtrs, ZGBTRS) (F77_CONST_CHAR_ARG_DECL,
00077 const octave_idx_type&, const octave_idx_type&,
00078 const octave_idx_type&, const octave_idx_type&,
00079 const Complex*, const octave_idx_type&,
00080 const octave_idx_type*, Complex*,
00081 const octave_idx_type&, octave_idx_type&
00082 F77_CHAR_ARG_LEN_DECL);
00083
00084 F77_RET_T
00085 F77_FUNC (zgbcon, ZGBCON) (F77_CONST_CHAR_ARG_DECL,
00086 const octave_idx_type&, const octave_idx_type&,
00087 const octave_idx_type&, Complex*,
00088 const octave_idx_type&, const octave_idx_type*,
00089 const double&, double&, Complex*, double*,
00090 octave_idx_type&
00091 F77_CHAR_ARG_LEN_DECL);
00092
00093 F77_RET_T
00094 F77_FUNC (zpbtrf, ZPBTRF) (F77_CONST_CHAR_ARG_DECL,
00095 const octave_idx_type&, const octave_idx_type&,
00096 Complex*, const octave_idx_type&, octave_idx_type&
00097 F77_CHAR_ARG_LEN_DECL);
00098
00099 F77_RET_T
00100 F77_FUNC (zpbtrs, ZPBTRS) (F77_CONST_CHAR_ARG_DECL,
00101 const octave_idx_type&, const octave_idx_type&,
00102 const octave_idx_type&, Complex*,
00103 const octave_idx_type&, Complex*,
00104 const octave_idx_type&, octave_idx_type&
00105 F77_CHAR_ARG_LEN_DECL);
00106
00107 F77_RET_T
00108 F77_FUNC (zpbcon, ZPBCON) (F77_CONST_CHAR_ARG_DECL,
00109 const octave_idx_type&, const octave_idx_type&,
00110 Complex*, const octave_idx_type&, const double&,
00111 double&, Complex*, double*, octave_idx_type&
00112 F77_CHAR_ARG_LEN_DECL);
00113
00114 F77_RET_T
00115 F77_FUNC (zgttrf, ZGTTRF) (const octave_idx_type&, Complex*, Complex*,
00116 Complex*, Complex*, octave_idx_type*,
00117 octave_idx_type&);
00118
00119 F77_RET_T
00120 F77_FUNC (zgttrs, ZGTTRS) (F77_CONST_CHAR_ARG_DECL,
00121 const octave_idx_type&, const octave_idx_type&,
00122 const Complex*, const Complex*, const Complex*,
00123 const Complex*, const octave_idx_type*,
00124 Complex *, const octave_idx_type&,
00125 octave_idx_type&
00126 F77_CHAR_ARG_LEN_DECL);
00127
00128 F77_RET_T
00129 F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&,
00130 double*, Complex*, Complex*,
00131 const octave_idx_type&, octave_idx_type&);
00132
00133 F77_RET_T
00134 F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&,
00135 Complex*, Complex*, Complex*, Complex*,
00136 const octave_idx_type&, octave_idx_type&);
00137 }
00138
00139 SparseComplexMatrix::SparseComplexMatrix (const SparseMatrix& a)
00140 : MSparse<Complex> (a)
00141 {
00142 }
00143
00144 SparseComplexMatrix::SparseComplexMatrix (const SparseBoolMatrix& a)
00145 : MSparse<Complex> (a.rows (), a.cols (), a.nnz ())
00146 {
00147 octave_idx_type nc = cols ();
00148 octave_idx_type nz = a.nnz ();
00149
00150 for (octave_idx_type i = 0; i < nc + 1; i++)
00151 cidx (i) = a.cidx (i);
00152
00153 for (octave_idx_type i = 0; i < nz; i++)
00154 {
00155 data (i) = Complex (a.data (i));
00156 ridx (i) = a.ridx (i);
00157 }
00158 }
00159
00160 SparseComplexMatrix::SparseComplexMatrix (const ComplexDiagMatrix& a)
00161 : MSparse<Complex> (a.rows (), a.cols (), a.length ())
00162 {
00163 octave_idx_type j = 0, l = a.length ();
00164 for (octave_idx_type i = 0; i < l; i++)
00165 {
00166 cidx (i) = j;
00167 if (a(i, i) != 0.0)
00168 {
00169 data (j) = a(i, i);
00170 ridx (j) = i;
00171 j++;
00172 }
00173 }
00174 for (octave_idx_type i = l; i <= a.cols (); i++)
00175 cidx(i) = j;
00176 }
00177 bool
00178 SparseComplexMatrix::operator == (const SparseComplexMatrix& a) const
00179 {
00180 octave_idx_type nr = rows ();
00181 octave_idx_type nc = cols ();
00182 octave_idx_type nz = nnz ();
00183 octave_idx_type nr_a = a.rows ();
00184 octave_idx_type nc_a = a.cols ();
00185 octave_idx_type nz_a = a.nnz ();
00186
00187 if (nr != nr_a || nc != nc_a || nz != nz_a)
00188 return false;
00189
00190 for (octave_idx_type i = 0; i < nc + 1; i++)
00191 if (cidx(i) != a.cidx(i))
00192 return false;
00193
00194 for (octave_idx_type i = 0; i < nz; i++)
00195 if (data(i) != a.data(i) || ridx(i) != a.ridx(i))
00196 return false;
00197
00198 return true;
00199 }
00200
00201 bool
00202 SparseComplexMatrix::operator != (const SparseComplexMatrix& a) const
00203 {
00204 return !(*this == a);
00205 }
00206
00207 bool
00208 SparseComplexMatrix::is_hermitian (void) const
00209 {
00210 octave_idx_type nr = rows ();
00211 octave_idx_type nc = cols ();
00212
00213 if (nr == nc && nr > 0)
00214 {
00215 for (octave_idx_type j = 0; j < nc; j++)
00216 {
00217 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00218 {
00219 octave_idx_type ri = ridx(i);
00220
00221 if (ri != j)
00222 {
00223 bool found = false;
00224
00225 for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++)
00226 {
00227 if (ridx(k) == j)
00228 {
00229 if (data(i) == conj(data(k)))
00230 found = true;
00231 break;
00232 }
00233 }
00234
00235 if (! found)
00236 return false;
00237 }
00238 }
00239 }
00240
00241 return true;
00242 }
00243
00244 return false;
00245 }
00246
00247 static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
00248
00249 SparseComplexMatrix
00250 SparseComplexMatrix::max (int dim) const
00251 {
00252 Array<octave_idx_type> dummy_idx;
00253 return max (dummy_idx, dim);
00254 }
00255
00256 SparseComplexMatrix
00257 SparseComplexMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
00258 {
00259 SparseComplexMatrix result;
00260 dim_vector dv = dims ();
00261
00262 if (dv.numel () == 0 || dim >= dv.length ())
00263 return result;
00264
00265 if (dim < 0)
00266 dim = dv.first_non_singleton ();
00267
00268 octave_idx_type nr = dv(0);
00269 octave_idx_type nc = dv(1);
00270
00271 if (dim == 0)
00272 {
00273 idx_arg.clear (1, nc);
00274 octave_idx_type nel = 0;
00275 for (octave_idx_type j = 0; j < nc; j++)
00276 {
00277 Complex tmp_max;
00278 double abs_max = octave_NaN;
00279 octave_idx_type idx_j = 0;
00280 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00281 {
00282 if (ridx(i) != idx_j)
00283 break;
00284 else
00285 idx_j++;
00286 }
00287
00288 if (idx_j != nr)
00289 {
00290 tmp_max = 0.;
00291 abs_max = 0.;
00292 }
00293
00294 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00295 {
00296 Complex tmp = data (i);
00297
00298 if (xisnan (tmp))
00299 continue;
00300
00301 double abs_tmp = std::abs (tmp);
00302
00303 if (xisnan (abs_max) || abs_tmp > abs_max)
00304 {
00305 idx_j = ridx (i);
00306 tmp_max = tmp;
00307 abs_max = abs_tmp;
00308 }
00309 }
00310
00311 idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_j;
00312 if (abs_max != 0.)
00313 nel++;
00314 }
00315
00316 result = SparseComplexMatrix (1, nc, nel);
00317
00318 octave_idx_type ii = 0;
00319 result.xcidx (0) = 0;
00320 for (octave_idx_type j = 0; j < nc; j++)
00321 {
00322 Complex tmp = elem (idx_arg(j), j);
00323 if (tmp != 0.)
00324 {
00325 result.xdata (ii) = tmp;
00326 result.xridx (ii++) = 0;
00327 }
00328 result.xcidx (j+1) = ii;
00329 }
00330 }
00331 else
00332 {
00333 idx_arg.resize (dim_vector (nr, 1), 0);
00334
00335 for (octave_idx_type i = cidx(0); i < cidx(1); i++)
00336 idx_arg.elem(ridx(i)) = -1;
00337
00338 for (octave_idx_type j = 0; j < nc; j++)
00339 for (octave_idx_type i = 0; i < nr; i++)
00340 {
00341 if (idx_arg.elem(i) != -1)
00342 continue;
00343 bool found = false;
00344 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
00345 if (ridx(k) == i)
00346 {
00347 found = true;
00348 break;
00349 }
00350
00351 if (!found)
00352 idx_arg.elem(i) = j;
00353
00354 }
00355
00356 for (octave_idx_type j = 0; j < nc; j++)
00357 {
00358 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00359 {
00360 octave_idx_type ir = ridx (i);
00361 octave_idx_type ix = idx_arg.elem (ir);
00362 Complex tmp = data (i);
00363
00364 if (xisnan (tmp))
00365 continue;
00366 else if (ix == -1 || std::abs(tmp) > std::abs(elem (ir, ix)))
00367 idx_arg.elem (ir) = j;
00368 }
00369 }
00370
00371 octave_idx_type nel = 0;
00372 for (octave_idx_type j = 0; j < nr; j++)
00373 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
00374 nel++;
00375
00376 result = SparseComplexMatrix (nr, 1, nel);
00377
00378 octave_idx_type ii = 0;
00379 result.xcidx (0) = 0;
00380 result.xcidx (1) = nel;
00381 for (octave_idx_type j = 0; j < nr; j++)
00382 {
00383 if (idx_arg(j) == -1)
00384 {
00385 idx_arg(j) = 0;
00386 result.xdata (ii) = Complex_NaN_result;
00387 result.xridx (ii++) = j;
00388 }
00389 else
00390 {
00391 Complex tmp = elem (j, idx_arg(j));
00392 if (tmp != 0.)
00393 {
00394 result.xdata (ii) = tmp;
00395 result.xridx (ii++) = j;
00396 }
00397 }
00398 }
00399 }
00400
00401 return result;
00402 }
00403
00404 SparseComplexMatrix
00405 SparseComplexMatrix::min (int dim) const
00406 {
00407 Array<octave_idx_type> dummy_idx;
00408 return min (dummy_idx, dim);
00409 }
00410
00411 SparseComplexMatrix
00412 SparseComplexMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
00413 {
00414 SparseComplexMatrix result;
00415 dim_vector dv = dims ();
00416
00417 if (dv.numel () == 0 || dim >= dv.length ())
00418 return result;
00419
00420 if (dim < 0)
00421 dim = dv.first_non_singleton ();
00422
00423 octave_idx_type nr = dv(0);
00424 octave_idx_type nc = dv(1);
00425
00426 if (dim == 0)
00427 {
00428 idx_arg.clear (1, nc);
00429 octave_idx_type nel = 0;
00430 for (octave_idx_type j = 0; j < nc; j++)
00431 {
00432 Complex tmp_min;
00433 double abs_min = octave_NaN;
00434 octave_idx_type idx_j = 0;
00435 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00436 {
00437 if (ridx(i) != idx_j)
00438 break;
00439 else
00440 idx_j++;
00441 }
00442
00443 if (idx_j != nr)
00444 {
00445 tmp_min = 0.;
00446 abs_min = 0.;
00447 }
00448
00449 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00450 {
00451 Complex tmp = data (i);
00452
00453 if (xisnan (tmp))
00454 continue;
00455
00456 double abs_tmp = std::abs (tmp);
00457
00458 if (xisnan (abs_min) || abs_tmp < abs_min)
00459 {
00460 idx_j = ridx (i);
00461 tmp_min = tmp;
00462 abs_min = abs_tmp;
00463 }
00464 }
00465
00466 idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_j;
00467 if (abs_min != 0.)
00468 nel++;
00469 }
00470
00471 result = SparseComplexMatrix (1, nc, nel);
00472
00473 octave_idx_type ii = 0;
00474 result.xcidx (0) = 0;
00475 for (octave_idx_type j = 0; j < nc; j++)
00476 {
00477 Complex tmp = elem (idx_arg(j), j);
00478 if (tmp != 0.)
00479 {
00480 result.xdata (ii) = tmp;
00481 result.xridx (ii++) = 0;
00482 }
00483 result.xcidx (j+1) = ii;
00484 }
00485 }
00486 else
00487 {
00488 idx_arg.resize (dim_vector (nr, 1), 0);
00489
00490 for (octave_idx_type i = cidx(0); i < cidx(1); i++)
00491 idx_arg.elem(ridx(i)) = -1;
00492
00493 for (octave_idx_type j = 0; j < nc; j++)
00494 for (octave_idx_type i = 0; i < nr; i++)
00495 {
00496 if (idx_arg.elem(i) != -1)
00497 continue;
00498 bool found = false;
00499 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
00500 if (ridx(k) == i)
00501 {
00502 found = true;
00503 break;
00504 }
00505
00506 if (!found)
00507 idx_arg.elem(i) = j;
00508
00509 }
00510
00511 for (octave_idx_type j = 0; j < nc; j++)
00512 {
00513 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00514 {
00515 octave_idx_type ir = ridx (i);
00516 octave_idx_type ix = idx_arg.elem (ir);
00517 Complex tmp = data (i);
00518
00519 if (xisnan (tmp))
00520 continue;
00521 else if (ix == -1 || std::abs(tmp) < std::abs(elem (ir, ix)))
00522 idx_arg.elem (ir) = j;
00523 }
00524 }
00525
00526 octave_idx_type nel = 0;
00527 for (octave_idx_type j = 0; j < nr; j++)
00528 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
00529 nel++;
00530
00531 result = SparseComplexMatrix (nr, 1, nel);
00532
00533 octave_idx_type ii = 0;
00534 result.xcidx (0) = 0;
00535 result.xcidx (1) = nel;
00536 for (octave_idx_type j = 0; j < nr; j++)
00537 {
00538 if (idx_arg(j) == -1)
00539 {
00540 idx_arg(j) = 0;
00541 result.xdata (ii) = Complex_NaN_result;
00542 result.xridx (ii++) = j;
00543 }
00544 else
00545 {
00546 Complex tmp = elem (j, idx_arg(j));
00547 if (tmp != 0.)
00548 {
00549 result.xdata (ii) = tmp;
00550 result.xridx (ii++) = j;
00551 }
00552 }
00553 }
00554 }
00555
00556 return result;
00557 }
00558
00559 ComplexRowVector
00560 SparseComplexMatrix::row (octave_idx_type i) const
00561 {
00562 octave_idx_type nc = columns ();
00563 ComplexRowVector retval (nc, 0);
00564
00565 for (octave_idx_type j = 0; j < nc; j++)
00566 for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
00567 {
00568 if (ridx (k) == i)
00569 {
00570 retval(j) = data (k);
00571 break;
00572 }
00573 }
00574
00575 return retval;
00576 }
00577
00578 ComplexColumnVector
00579 SparseComplexMatrix::column (octave_idx_type i) const
00580 {
00581 octave_idx_type nr = rows ();
00582 ComplexColumnVector retval (nr, 0);
00583
00584 for (octave_idx_type k = cidx (i); k < cidx (i+1); k++)
00585 retval(ridx (k)) = data (k);
00586
00587 return retval;
00588 }
00589
00590
00591
00592 SparseComplexMatrix&
00593 SparseComplexMatrix::insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c)
00594 {
00595 SparseComplexMatrix tmp (a);
00596 return insert (tmp , r, c);
00597 }
00598
00599 SparseComplexMatrix&
00600 SparseComplexMatrix::insert (const SparseComplexMatrix& a, octave_idx_type r, octave_idx_type c)
00601 {
00602 MSparse<Complex>::insert (a, r, c);
00603 return *this;
00604 }
00605
00606 SparseComplexMatrix&
00607 SparseComplexMatrix::insert (const SparseMatrix& a, const Array<octave_idx_type>& indx)
00608 {
00609 SparseComplexMatrix tmp (a);
00610 return insert (tmp , indx);
00611 }
00612
00613 SparseComplexMatrix&
00614 SparseComplexMatrix::insert (const SparseComplexMatrix& a, const Array<octave_idx_type>& indx)
00615 {
00616 MSparse<Complex>::insert (a, indx);
00617 return *this;
00618 }
00619
00620 SparseComplexMatrix
00621 SparseComplexMatrix::concat (const SparseComplexMatrix& rb,
00622 const Array<octave_idx_type>& ra_idx)
00623 {
00624
00625 if (rb.rows () > 0 && rb.cols () > 0)
00626 insert (rb, ra_idx(0), ra_idx(1));
00627 return *this;
00628 }
00629
00630 SparseComplexMatrix
00631 SparseComplexMatrix::concat (const SparseMatrix& rb, const Array<octave_idx_type>& ra_idx)
00632 {
00633 SparseComplexMatrix tmp (rb);
00634 if (rb.rows () > 0 && rb.cols () > 0)
00635 insert (tmp, ra_idx(0), ra_idx(1));
00636 return *this;
00637 }
00638
00639 ComplexMatrix
00640 SparseComplexMatrix::matrix_value (void) const
00641 {
00642 return Sparse<Complex>::array_value ();
00643 }
00644
00645 SparseComplexMatrix
00646 SparseComplexMatrix::hermitian (void) const
00647 {
00648 octave_idx_type nr = rows ();
00649 octave_idx_type nc = cols ();
00650 octave_idx_type nz = nnz ();
00651 SparseComplexMatrix retval (nc, nr, nz);
00652
00653 for (octave_idx_type i = 0; i < nz; i++)
00654 retval.xcidx (ridx (i) + 1)++;
00655
00656 nz = 0;
00657 for (octave_idx_type i = 1; i <= nr; i++)
00658 {
00659 const octave_idx_type tmp = retval.xcidx (i);
00660 retval.xcidx (i) = nz;
00661 nz += tmp;
00662 }
00663
00664
00665 for (octave_idx_type j = 0; j < nc; j++)
00666 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
00667 {
00668 octave_idx_type q = retval.xcidx (ridx (k) + 1)++;
00669 retval.xridx (q) = j;
00670 retval.xdata (q) = conj (data (k));
00671 }
00672 assert (nnz () == retval.xcidx (nr));
00673
00674
00675
00676 return retval;
00677 }
00678
00679 SparseComplexMatrix
00680 conj (const SparseComplexMatrix& a)
00681 {
00682 octave_idx_type nr = a.rows ();
00683 octave_idx_type nc = a.cols ();
00684 octave_idx_type nz = a.nnz ();
00685 SparseComplexMatrix retval (nc, nr, nz);
00686
00687 for (octave_idx_type i = 0; i < nc + 1; i++)
00688 retval.cidx (i) = a.cidx (i);
00689
00690 for (octave_idx_type i = 0; i < nz; i++)
00691 {
00692 retval.data (i) = conj (a.data (i));
00693 retval.ridx (i) = a.ridx (i);
00694 }
00695
00696 return retval;
00697 }
00698
00699 SparseComplexMatrix
00700 SparseComplexMatrix::inverse (void) const
00701 {
00702 octave_idx_type info;
00703 double rcond;
00704 MatrixType mattype (*this);
00705 return inverse (mattype, info, rcond, 0, 0);
00706 }
00707
00708 SparseComplexMatrix
00709 SparseComplexMatrix::inverse (MatrixType& mattype) const
00710 {
00711 octave_idx_type info;
00712 double rcond;
00713 return inverse (mattype, info, rcond, 0, 0);
00714 }
00715
00716 SparseComplexMatrix
00717 SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info) const
00718 {
00719 double rcond;
00720 return inverse (mattype, info, rcond, 0, 0);
00721 }
00722
00723 SparseComplexMatrix
00724 SparseComplexMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info,
00725 double& rcond, const bool,
00726 const bool calccond) const
00727 {
00728 SparseComplexMatrix retval;
00729
00730 octave_idx_type nr = rows ();
00731 octave_idx_type nc = cols ();
00732 info = 0;
00733
00734 if (nr == 0 || nc == 0 || nr != nc)
00735 (*current_liboctave_error_handler) ("inverse requires square matrix");
00736 else
00737 {
00738
00739 int typ = mattyp.type ();
00740 mattyp.info ();
00741
00742 if (typ == MatrixType::Diagonal ||
00743 typ == MatrixType::Permuted_Diagonal)
00744 {
00745 if (typ == MatrixType::Permuted_Diagonal)
00746 retval = transpose();
00747 else
00748 retval = *this;
00749
00750
00751 Complex *v = retval.data();
00752
00753 if (calccond)
00754 {
00755 double dmax = 0., dmin = octave_Inf;
00756 for (octave_idx_type i = 0; i < nr; i++)
00757 {
00758 double tmp = std::abs(v[i]);
00759 if (tmp > dmax)
00760 dmax = tmp;
00761 if (tmp < dmin)
00762 dmin = tmp;
00763 }
00764 rcond = dmin / dmax;
00765 }
00766
00767 for (octave_idx_type i = 0; i < nr; i++)
00768 v[i] = 1.0 / v[i];
00769 }
00770 else
00771 (*current_liboctave_error_handler) ("incorrect matrix type");
00772 }
00773
00774 return retval;
00775 }
00776
00777 SparseComplexMatrix
00778 SparseComplexMatrix::tinverse (MatrixType &mattyp, octave_idx_type& info,
00779 double& rcond, const bool,
00780 const bool calccond) const
00781 {
00782 SparseComplexMatrix retval;
00783
00784 octave_idx_type nr = rows ();
00785 octave_idx_type nc = cols ();
00786 info = 0;
00787
00788 if (nr == 0 || nc == 0 || nr != nc)
00789 (*current_liboctave_error_handler) ("inverse requires square matrix");
00790 else
00791 {
00792
00793 int typ = mattyp.type ();
00794 mattyp.info ();
00795
00796 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper ||
00797 typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
00798 {
00799 double anorm = 0.;
00800 double ainvnorm = 0.;
00801
00802 if (calccond)
00803 {
00804
00805 for (octave_idx_type j = 0; j < nr; j++)
00806 {
00807 double atmp = 0.;
00808 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
00809 atmp += std::abs(data(i));
00810 if (atmp > anorm)
00811 anorm = atmp;
00812 }
00813 }
00814
00815 if (typ == MatrixType::Upper || typ == MatrixType::Lower)
00816 {
00817 octave_idx_type nz = nnz ();
00818 octave_idx_type cx = 0;
00819 octave_idx_type nz2 = nz;
00820 retval = SparseComplexMatrix (nr, nc, nz2);
00821
00822 for (octave_idx_type i = 0; i < nr; i++)
00823 {
00824 octave_quit ();
00825
00826 octave_idx_type cx_colstart = cx;
00827
00828 if (cx == nz2)
00829 {
00830 nz2 *= 2;
00831 retval.change_capacity (nz2);
00832 }
00833
00834 retval.xcidx(i) = cx;
00835 retval.xridx(cx) = i;
00836 retval.xdata(cx) = 1.0;
00837 cx++;
00838
00839
00840 for (octave_idx_type j = i+1; j < nr; j++)
00841 {
00842 Complex v = 0.;
00843
00844 octave_idx_type colXp = retval.xcidx(i);
00845 octave_idx_type colUp = cidx(j);
00846 octave_idx_type rpX, rpU;
00847
00848 if (cidx(j) == cidx(j+1))
00849 {
00850 (*current_liboctave_error_handler)
00851 ("division by zero");
00852 goto inverse_singular;
00853 }
00854
00855 do
00856 {
00857 octave_quit ();
00858 rpX = retval.xridx(colXp);
00859 rpU = ridx(colUp);
00860
00861 if (rpX < rpU)
00862 colXp++;
00863 else if (rpX > rpU)
00864 colUp++;
00865 else
00866 {
00867 v -= retval.xdata(colXp) * data(colUp);
00868 colXp++;
00869 colUp++;
00870 }
00871 } while ((rpX<j) && (rpU<j) &&
00872 (colXp<cx) && (colUp<nz));
00873
00874
00875
00876 if (typ == MatrixType::Upper)
00877 colUp = cidx(j+1) - 1;
00878 else
00879 colUp = cidx(j);
00880 Complex pivot = data(colUp);
00881 if (pivot == 0. || ridx(colUp) != j)
00882 {
00883 (*current_liboctave_error_handler)
00884 ("division by zero");
00885 goto inverse_singular;
00886 }
00887
00888 if (v != 0.)
00889 {
00890 if (cx == nz2)
00891 {
00892 nz2 *= 2;
00893 retval.change_capacity (nz2);
00894 }
00895
00896 retval.xridx(cx) = j;
00897 retval.xdata(cx) = v / pivot;
00898 cx++;
00899 }
00900 }
00901
00902
00903 octave_idx_type colUp;
00904 if (typ == MatrixType::Upper)
00905 colUp = cidx(i+1) - 1;
00906 else
00907 colUp = cidx(i);
00908 Complex pivot = data(colUp);
00909 if (pivot == 0. || ridx(colUp) != i)
00910 {
00911 (*current_liboctave_error_handler) ("division by zero");
00912 goto inverse_singular;
00913 }
00914
00915 if (pivot != 1.0)
00916 for (octave_idx_type j = cx_colstart; j < cx; j++)
00917 retval.xdata(j) /= pivot;
00918 }
00919 retval.xcidx(nr) = cx;
00920 retval.maybe_compress ();
00921 }
00922 else
00923 {
00924 octave_idx_type nz = nnz ();
00925 octave_idx_type cx = 0;
00926 octave_idx_type nz2 = nz;
00927 retval = SparseComplexMatrix (nr, nc, nz2);
00928
00929 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
00930 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
00931
00932 octave_idx_type *perm = mattyp.triangular_perm();
00933 if (typ == MatrixType::Permuted_Upper)
00934 {
00935 for (octave_idx_type i = 0; i < nr; i++)
00936 rperm[perm[i]] = i;
00937 }
00938 else
00939 {
00940 for (octave_idx_type i = 0; i < nr; i++)
00941 rperm[i] = perm[i];
00942 for (octave_idx_type i = 0; i < nr; i++)
00943 perm[rperm[i]] = i;
00944 }
00945
00946 for (octave_idx_type i = 0; i < nr; i++)
00947 {
00948 octave_quit ();
00949 octave_idx_type iidx = rperm[i];
00950
00951 for (octave_idx_type j = 0; j < nr; j++)
00952 work[j] = 0.;
00953
00954
00955 work[iidx] = 1.0;
00956
00957
00958 for (octave_idx_type j = iidx+1; j < nr; j++)
00959 {
00960 Complex v = 0.;
00961 octave_idx_type jidx = perm[j];
00962
00963 for (octave_idx_type k = cidx(jidx);
00964 k < cidx(jidx+1); k++)
00965 {
00966 octave_quit ();
00967 v -= work[ridx(k)] * data(k);
00968 }
00969
00970
00971 Complex pivot;
00972 if (typ == MatrixType::Permuted_Upper)
00973 pivot = data(cidx(jidx+1) - 1);
00974 else
00975 pivot = data(cidx(jidx));
00976 if (pivot == 0.)
00977 {
00978 (*current_liboctave_error_handler)
00979 ("division by zero");
00980 goto inverse_singular;
00981 }
00982
00983 work[j] = v / pivot;
00984 }
00985
00986
00987 octave_idx_type colUp;
00988 if (typ == MatrixType::Permuted_Upper)
00989 colUp = cidx(perm[iidx]+1) - 1;
00990 else
00991 colUp = cidx(perm[iidx]);
00992
00993 Complex pivot = data(colUp);
00994 if (pivot == 0.)
00995 {
00996 (*current_liboctave_error_handler)
00997 ("division by zero");
00998 goto inverse_singular;
00999 }
01000
01001 octave_idx_type new_cx = cx;
01002 for (octave_idx_type j = iidx; j < nr; j++)
01003 if (work[j] != 0.0)
01004 {
01005 new_cx++;
01006 if (pivot != 1.0)
01007 work[j] /= pivot;
01008 }
01009
01010 if (cx < new_cx)
01011 {
01012 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
01013 retval.change_capacity (nz2);
01014 }
01015
01016 retval.xcidx(i) = cx;
01017 for (octave_idx_type j = iidx; j < nr; j++)
01018 if (work[j] != 0.)
01019 {
01020 retval.xridx(cx) = j;
01021 retval.xdata(cx++) = work[j];
01022 }
01023 }
01024
01025 retval.xcidx(nr) = cx;
01026 retval.maybe_compress ();
01027 }
01028
01029 if (calccond)
01030 {
01031
01032 for (octave_idx_type j = 0; j < nr; j++)
01033 {
01034 double atmp = 0.;
01035 for (octave_idx_type i = retval.cidx(j);
01036 i < retval.cidx(j+1); i++)
01037 atmp += std::abs(retval.data(i));
01038 if (atmp > ainvnorm)
01039 ainvnorm = atmp;
01040 }
01041
01042 rcond = 1. / ainvnorm / anorm;
01043 }
01044 }
01045 else
01046 (*current_liboctave_error_handler) ("incorrect matrix type");
01047 }
01048
01049 return retval;
01050
01051 inverse_singular:
01052 return SparseComplexMatrix();
01053 }
01054
01055 SparseComplexMatrix
01056 SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info,
01057 double& rcond, int, int calc_cond) const
01058 {
01059 int typ = mattype.type (false);
01060 SparseComplexMatrix ret;
01061
01062 if (typ == MatrixType::Unknown)
01063 typ = mattype.type (*this);
01064
01065 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
01066 ret = dinverse (mattype, info, rcond, true, calc_cond);
01067 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
01068 ret = tinverse (mattype, info, rcond, true, calc_cond).transpose();
01069 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
01070 {
01071 MatrixType newtype = mattype.transpose();
01072 ret = transpose().tinverse (newtype, info, rcond, true, calc_cond);
01073 }
01074 else
01075 {
01076 if (mattype.is_hermitian())
01077 {
01078 MatrixType tmp_typ (MatrixType::Upper);
01079 SparseComplexCHOL fact (*this, info, false);
01080 rcond = fact.rcond();
01081 if (info == 0)
01082 {
01083 double rcond2;
01084 SparseMatrix Q = fact.Q();
01085 SparseComplexMatrix InvL = fact.L().transpose().
01086 tinverse(tmp_typ, info, rcond2, true, false);
01087 ret = Q * InvL.hermitian() * InvL * Q.transpose();
01088 }
01089 else
01090 {
01091
01092 mattype.mark_as_unsymmetric ();
01093 typ = MatrixType::Full;
01094 }
01095 }
01096
01097 if (!mattype.is_hermitian())
01098 {
01099 octave_idx_type n = rows();
01100 ColumnVector Qinit(n);
01101 for (octave_idx_type i = 0; i < n; i++)
01102 Qinit(i) = i;
01103
01104 MatrixType tmp_typ (MatrixType::Upper);
01105 SparseComplexLU fact (*this, Qinit, Matrix (), false, false);
01106 rcond = fact.rcond();
01107 double rcond2;
01108 SparseComplexMatrix InvL = fact.L().transpose().
01109 tinverse(tmp_typ, info, rcond2, true, false);
01110 SparseComplexMatrix InvU = fact.U().
01111 tinverse(tmp_typ, info, rcond2, true, false).transpose();
01112 ret = fact.Pc().transpose() * InvU * InvL * fact.Pr();
01113 }
01114 }
01115
01116 return ret;
01117 }
01118
01119 ComplexDET
01120 SparseComplexMatrix::determinant (void) const
01121 {
01122 octave_idx_type info;
01123 double rcond;
01124 return determinant (info, rcond, 0);
01125 }
01126
01127 ComplexDET
01128 SparseComplexMatrix::determinant (octave_idx_type& info) const
01129 {
01130 double rcond;
01131 return determinant (info, rcond, 0);
01132 }
01133
01134 ComplexDET
01135 SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, int) const
01136 {
01137 ComplexDET retval;
01138 #ifdef HAVE_UMFPACK
01139
01140 octave_idx_type nr = rows ();
01141 octave_idx_type nc = cols ();
01142
01143 if (nr == 0 || nc == 0 || nr != nc)
01144 {
01145 retval = ComplexDET (1.0);
01146 }
01147 else
01148 {
01149 err = 0;
01150
01151
01152 Matrix Control (UMFPACK_CONTROL, 1);
01153 double *control = Control.fortran_vec ();
01154 UMFPACK_ZNAME (defaults) (control);
01155
01156 double tmp = octave_sparse_params::get_key ("spumoni");
01157 if (!xisnan (tmp))
01158 Control (UMFPACK_PRL) = tmp;
01159
01160 tmp = octave_sparse_params::get_key ("piv_tol");
01161 if (!xisnan (tmp))
01162 {
01163 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
01164 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
01165 }
01166
01167
01168 tmp = octave_sparse_params::get_key ("autoamd");
01169 if (!xisnan (tmp))
01170 Control (UMFPACK_FIXQ) = tmp;
01171
01172
01173 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
01174
01175 UMFPACK_ZNAME (report_control) (control);
01176
01177 const octave_idx_type *Ap = cidx ();
01178 const octave_idx_type *Ai = ridx ();
01179 const Complex *Ax = data ();
01180
01181 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
01182 reinterpret_cast<const double *> (Ax),
01183 0, 1, control);
01184
01185 void *Symbolic;
01186 Matrix Info (1, UMFPACK_INFO);
01187 double *info = Info.fortran_vec ();
01188 int status = UMFPACK_ZNAME (qsymbolic)
01189 (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0,
01190 0, &Symbolic, control, info);
01191
01192 if (status < 0)
01193 {
01194 (*current_liboctave_error_handler)
01195 ("SparseComplexMatrix::determinant symbolic factorization failed");
01196
01197 UMFPACK_ZNAME (report_status) (control, status);
01198 UMFPACK_ZNAME (report_info) (control, info);
01199
01200 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
01201 }
01202 else
01203 {
01204 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
01205
01206 void *Numeric;
01207 status
01208 = UMFPACK_ZNAME (numeric) (Ap, Ai,
01209 reinterpret_cast<const double *> (Ax),
01210 0, Symbolic, &Numeric, control, info) ;
01211 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
01212
01213 rcond = Info (UMFPACK_RCOND);
01214
01215 if (status < 0)
01216 {
01217 (*current_liboctave_error_handler)
01218 ("SparseComplexMatrix::determinant numeric factorization failed");
01219
01220 UMFPACK_ZNAME (report_status) (control, status);
01221 UMFPACK_ZNAME (report_info) (control, info);
01222
01223 UMFPACK_ZNAME (free_numeric) (&Numeric);
01224 }
01225 else
01226 {
01227 UMFPACK_ZNAME (report_numeric) (Numeric, control);
01228
01229 double c10[2], e10;
01230
01231 status = UMFPACK_ZNAME (get_determinant) (c10, 0, &e10,
01232 Numeric, info);
01233
01234 if (status < 0)
01235 {
01236 (*current_liboctave_error_handler)
01237 ("SparseComplexMatrix::determinant error calculating determinant");
01238
01239 UMFPACK_ZNAME (report_status) (control, status);
01240 UMFPACK_ZNAME (report_info) (control, info);
01241 }
01242 else
01243 retval = ComplexDET (Complex (c10[0], c10[1]), e10, 10);
01244
01245 UMFPACK_ZNAME (free_numeric) (&Numeric);
01246 }
01247 }
01248 }
01249 #else
01250 (*current_liboctave_error_handler) ("UMFPACK not installed");
01251 #endif
01252
01253 return retval;
01254 }
01255
01256 ComplexMatrix
01257 SparseComplexMatrix::dsolve (MatrixType &mattype, const Matrix& b,
01258 octave_idx_type& err, double& rcond,
01259 solve_singularity_handler, bool calc_cond) const
01260 {
01261 ComplexMatrix retval;
01262
01263 octave_idx_type nr = rows ();
01264 octave_idx_type nc = cols ();
01265 octave_idx_type nm = (nc < nr ? nc : nr);
01266 err = 0;
01267
01268 if (nr != b.rows ())
01269 (*current_liboctave_error_handler)
01270 ("matrix dimension mismatch solution of linear equations");
01271 else if (nr == 0 || nc == 0 || b.cols () == 0)
01272 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
01273 else
01274 {
01275
01276 int typ = mattype.type ();
01277 mattype.info ();
01278
01279 if (typ == MatrixType::Diagonal ||
01280 typ == MatrixType::Permuted_Diagonal)
01281 {
01282 retval.resize (nc, b.cols(), Complex(0.,0.));
01283 if (typ == MatrixType::Diagonal)
01284 for (octave_idx_type j = 0; j < b.cols(); j++)
01285 for (octave_idx_type i = 0; i < nm; i++)
01286 retval(i,j) = b(i,j) / data (i);
01287 else
01288 for (octave_idx_type j = 0; j < b.cols(); j++)
01289 for (octave_idx_type k = 0; k < nc; k++)
01290 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
01291 retval(k,j) = b(ridx(i),j) / data (i);
01292
01293 if (calc_cond)
01294 {
01295 double dmax = 0., dmin = octave_Inf;
01296 for (octave_idx_type i = 0; i < nm; i++)
01297 {
01298 double tmp = std::abs(data(i));
01299 if (tmp > dmax)
01300 dmax = tmp;
01301 if (tmp < dmin)
01302 dmin = tmp;
01303 }
01304 rcond = dmin / dmax;
01305 }
01306 else
01307 rcond = 1.0;
01308 }
01309 else
01310 (*current_liboctave_error_handler) ("incorrect matrix type");
01311 }
01312
01313 return retval;
01314 }
01315
01316 SparseComplexMatrix
01317 SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseMatrix& b,
01318 octave_idx_type& err, double& rcond,
01319 solve_singularity_handler,
01320 bool calc_cond) const
01321 {
01322 SparseComplexMatrix retval;
01323
01324 octave_idx_type nr = rows ();
01325 octave_idx_type nc = cols ();
01326 octave_idx_type nm = (nc < nr ? nc : nr);
01327 err = 0;
01328
01329 if (nr != b.rows ())
01330 (*current_liboctave_error_handler)
01331 ("matrix dimension mismatch solution of linear equations");
01332 else if (nr == 0 || nc == 0 || b.cols () == 0)
01333 retval = SparseComplexMatrix (nc, b.cols ());
01334 else
01335 {
01336
01337 int typ = mattype.type ();
01338 mattype.info ();
01339
01340 if (typ == MatrixType::Diagonal ||
01341 typ == MatrixType::Permuted_Diagonal)
01342 {
01343 octave_idx_type b_nc = b.cols ();
01344 octave_idx_type b_nz = b.nnz ();
01345 retval = SparseComplexMatrix (nc, b_nc, b_nz);
01346
01347 retval.xcidx(0) = 0;
01348 octave_idx_type ii = 0;
01349 if (typ == MatrixType::Diagonal)
01350 for (octave_idx_type j = 0; j < b.cols(); j++)
01351 {
01352 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
01353 {
01354 if (b.ridx(i) >= nm)
01355 break;
01356 retval.xridx (ii) = b.ridx(i);
01357 retval.xdata (ii++) = b.data(i) / data (b.ridx (i));
01358 }
01359 retval.xcidx(j+1) = ii;
01360 }
01361 else
01362 for (octave_idx_type j = 0; j < b.cols(); j++)
01363 {
01364 for (octave_idx_type l = 0; l < nc; l++)
01365 for (octave_idx_type i = cidx(l); i < cidx(l+1); i++)
01366 {
01367 bool found = false;
01368 octave_idx_type k;
01369 for (k = b.cidx(j); k < b.cidx(j+1); k++)
01370 if (ridx(i) == b.ridx(k))
01371 {
01372 found = true;
01373 break;
01374 }
01375 if (found)
01376 {
01377 retval.xridx (ii) = l;
01378 retval.xdata (ii++) = b.data(k) / data (i);
01379 }
01380 }
01381 retval.xcidx(j+1) = ii;
01382 }
01383
01384 if (calc_cond)
01385 {
01386 double dmax = 0., dmin = octave_Inf;
01387 for (octave_idx_type i = 0; i < nm; i++)
01388 {
01389 double tmp = std::abs(data(i));
01390 if (tmp > dmax)
01391 dmax = tmp;
01392 if (tmp < dmin)
01393 dmin = tmp;
01394 }
01395 rcond = dmin / dmax;
01396 }
01397 else
01398 rcond = 1.0;
01399 }
01400 else
01401 (*current_liboctave_error_handler) ("incorrect matrix type");
01402 }
01403
01404 return retval;
01405 }
01406
01407 ComplexMatrix
01408 SparseComplexMatrix::dsolve (MatrixType &mattype, const ComplexMatrix& b,
01409 octave_idx_type& err, double& rcond,
01410 solve_singularity_handler,
01411 bool calc_cond) const
01412 {
01413 ComplexMatrix retval;
01414
01415 octave_idx_type nr = rows ();
01416 octave_idx_type nc = cols ();
01417 octave_idx_type nm = (nc < nr ? nc : nr);
01418 err = 0;
01419
01420 if (nr != b.rows ())
01421 (*current_liboctave_error_handler)
01422 ("matrix dimension mismatch solution of linear equations");
01423 else if (nr == 0 || nc == 0 || b.cols () == 0)
01424 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
01425 else
01426 {
01427
01428 int typ = mattype.type ();
01429 mattype.info ();
01430
01431 if (typ == MatrixType::Diagonal ||
01432 typ == MatrixType::Permuted_Diagonal)
01433 {
01434 retval.resize (nc, b.cols(), Complex(0.,0.));
01435 if (typ == MatrixType::Diagonal)
01436 for (octave_idx_type j = 0; j < b.cols(); j++)
01437 for (octave_idx_type i = 0; i < nm; i++)
01438 retval(i,j) = b(i,j) / data (i);
01439 else
01440 for (octave_idx_type j = 0; j < b.cols(); j++)
01441 for (octave_idx_type k = 0; k < nc; k++)
01442 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
01443 retval(k,j) = b(ridx(i),j) / data (i);
01444
01445 if (calc_cond)
01446 {
01447 double dmax = 0., dmin = octave_Inf;
01448 for (octave_idx_type i = 0; i < nr; i++)
01449 {
01450 double tmp = std::abs(data(i));
01451 if (tmp > dmax)
01452 dmax = tmp;
01453 if (tmp < dmin)
01454 dmin = tmp;
01455 }
01456 rcond = dmin / dmax;
01457 }
01458 else
01459 rcond = 1.0;
01460 }
01461 else
01462 (*current_liboctave_error_handler) ("incorrect matrix type");
01463 }
01464
01465 return retval;
01466 }
01467
01468 SparseComplexMatrix
01469 SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseComplexMatrix& b,
01470 octave_idx_type& err, double& rcond,
01471 solve_singularity_handler,
01472 bool calc_cond) const
01473 {
01474 SparseComplexMatrix retval;
01475
01476 octave_idx_type nr = rows ();
01477 octave_idx_type nc = cols ();
01478 octave_idx_type nm = (nc < nr ? nc : nr);
01479 err = 0;
01480
01481 if (nr != b.rows ())
01482 (*current_liboctave_error_handler)
01483 ("matrix dimension mismatch solution of linear equations");
01484 else if (nr == 0 || nc == 0 || b.cols () == 0)
01485 retval = SparseComplexMatrix (nc, b.cols ());
01486 else
01487 {
01488
01489 int typ = mattype.type ();
01490 mattype.info ();
01491
01492 if (typ == MatrixType::Diagonal ||
01493 typ == MatrixType::Permuted_Diagonal)
01494 {
01495 octave_idx_type b_nc = b.cols ();
01496 octave_idx_type b_nz = b.nnz ();
01497 retval = SparseComplexMatrix (nc, b_nc, b_nz);
01498
01499 retval.xcidx(0) = 0;
01500 octave_idx_type ii = 0;
01501 if (typ == MatrixType::Diagonal)
01502 for (octave_idx_type j = 0; j < b.cols(); j++)
01503 {
01504 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
01505 {
01506 if (b.ridx(i) >= nm)
01507 break;
01508 retval.xridx (ii) = b.ridx(i);
01509 retval.xdata (ii++) = b.data(i) / data (b.ridx (i));
01510 }
01511 retval.xcidx(j+1) = ii;
01512 }
01513 else
01514 for (octave_idx_type j = 0; j < b.cols(); j++)
01515 {
01516 for (octave_idx_type l = 0; l < nc; l++)
01517 for (octave_idx_type i = cidx(l); i < cidx(l+1); i++)
01518 {
01519 bool found = false;
01520 octave_idx_type k;
01521 for (k = b.cidx(j); k < b.cidx(j+1); k++)
01522 if (ridx(i) == b.ridx(k))
01523 {
01524 found = true;
01525 break;
01526 }
01527 if (found)
01528 {
01529 retval.xridx (ii) = l;
01530 retval.xdata (ii++) = b.data(k) / data (i);
01531 }
01532 }
01533 retval.xcidx(j+1) = ii;
01534 }
01535
01536 if (calc_cond)
01537 {
01538 double dmax = 0., dmin = octave_Inf;
01539 for (octave_idx_type i = 0; i < nm; i++)
01540 {
01541 double tmp = std::abs(data(i));
01542 if (tmp > dmax)
01543 dmax = tmp;
01544 if (tmp < dmin)
01545 dmin = tmp;
01546 }
01547 rcond = dmin / dmax;
01548 }
01549 else
01550 rcond = 1.0;
01551 }
01552 else
01553 (*current_liboctave_error_handler) ("incorrect matrix type");
01554 }
01555
01556 return retval;
01557 }
01558
01559 ComplexMatrix
01560 SparseComplexMatrix::utsolve (MatrixType &mattype, const Matrix& b,
01561 octave_idx_type& err, double& rcond,
01562 solve_singularity_handler sing_handler,
01563 bool calc_cond) const
01564 {
01565 ComplexMatrix retval;
01566
01567 octave_idx_type nr = rows ();
01568 octave_idx_type nc = cols ();
01569 octave_idx_type nm = (nc > nr ? nc : nr);
01570 err = 0;
01571
01572 if (nr != b.rows ())
01573 (*current_liboctave_error_handler)
01574 ("matrix dimension mismatch solution of linear equations");
01575 else if (nr == 0 || nc == 0 || b.cols () == 0)
01576 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
01577 else
01578 {
01579
01580 int typ = mattype.type ();
01581 mattype.info ();
01582
01583 if (typ == MatrixType::Permuted_Upper ||
01584 typ == MatrixType::Upper)
01585 {
01586 double anorm = 0.;
01587 double ainvnorm = 0.;
01588 octave_idx_type b_nc = b.cols ();
01589 rcond = 1.;
01590
01591 if (calc_cond)
01592 {
01593
01594 for (octave_idx_type j = 0; j < nc; j++)
01595 {
01596 double atmp = 0.;
01597 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
01598 atmp += std::abs(data(i));
01599 if (atmp > anorm)
01600 anorm = atmp;
01601 }
01602 }
01603
01604 if (typ == MatrixType::Permuted_Upper)
01605 {
01606 retval.resize (nc, b_nc);
01607 octave_idx_type *perm = mattype.triangular_perm ();
01608 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
01609
01610 for (octave_idx_type j = 0; j < b_nc; j++)
01611 {
01612 for (octave_idx_type i = 0; i < nr; i++)
01613 work[i] = b(i,j);
01614 for (octave_idx_type i = nr; i < nc; i++)
01615 work[i] = 0.;
01616
01617 for (octave_idx_type k = nc-1; k >= 0; k--)
01618 {
01619 octave_idx_type kidx = perm[k];
01620
01621 if (work[k] != 0.)
01622 {
01623 if (ridx(cidx(kidx+1)-1) != k ||
01624 data(cidx(kidx+1)-1) == 0.)
01625 {
01626 err = -2;
01627 goto triangular_error;
01628 }
01629
01630 Complex tmp = work[k] / data(cidx(kidx+1)-1);
01631 work[k] = tmp;
01632 for (octave_idx_type i = cidx(kidx);
01633 i < cidx(kidx+1)-1; i++)
01634 {
01635 octave_idx_type iidx = ridx(i);
01636 work[iidx] = work[iidx] - tmp * data(i);
01637 }
01638 }
01639 }
01640
01641 for (octave_idx_type i = 0; i < nc; i++)
01642 retval (perm[i], j) = work[i];
01643 }
01644
01645 if (calc_cond)
01646 {
01647
01648 for (octave_idx_type i = 0; i < nm; i++)
01649 work[i] = 0.;
01650
01651 for (octave_idx_type j = 0; j < nr; j++)
01652 {
01653 work[j] = 1.;
01654
01655 for (octave_idx_type k = j; k >= 0; k--)
01656 {
01657 octave_idx_type iidx = perm[k];
01658
01659 if (work[k] != 0.)
01660 {
01661 Complex tmp = work[k] / data(cidx(iidx+1)-1);
01662 work[k] = tmp;
01663 for (octave_idx_type i = cidx(iidx);
01664 i < cidx(iidx+1)-1; i++)
01665 {
01666 octave_idx_type idx2 = ridx(i);
01667 work[idx2] = work[idx2] - tmp * data(i);
01668 }
01669 }
01670 }
01671 double atmp = 0;
01672 for (octave_idx_type i = 0; i < j+1; i++)
01673 {
01674 atmp += std::abs(work[i]);
01675 work[i] = 0.;
01676 }
01677 if (atmp > ainvnorm)
01678 ainvnorm = atmp;
01679 }
01680 rcond = 1. / ainvnorm / anorm;
01681 }
01682 }
01683 else
01684 {
01685 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
01686 retval.resize (nc, b_nc);
01687
01688 for (octave_idx_type j = 0; j < b_nc; j++)
01689 {
01690 for (octave_idx_type i = 0; i < nr; i++)
01691 work[i] = b(i,j);
01692 for (octave_idx_type i = nr; i < nc; i++)
01693 work[i] = 0.;
01694
01695 for (octave_idx_type k = nc-1; k >= 0; k--)
01696 {
01697 if (work[k] != 0.)
01698 {
01699 if (ridx(cidx(k+1)-1) != k ||
01700 data(cidx(k+1)-1) == 0.)
01701 {
01702 err = -2;
01703 goto triangular_error;
01704 }
01705
01706 Complex tmp = work[k] / data(cidx(k+1)-1);
01707 work[k] = tmp;
01708 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
01709 {
01710 octave_idx_type iidx = ridx(i);
01711 work[iidx] = work[iidx] - tmp * data(i);
01712 }
01713 }
01714 }
01715
01716 for (octave_idx_type i = 0; i < nc; i++)
01717 retval.xelem (i, j) = work[i];
01718 }
01719
01720 if (calc_cond)
01721 {
01722
01723 for (octave_idx_type i = 0; i < nm; i++)
01724 work[i] = 0.;
01725
01726 for (octave_idx_type j = 0; j < nr; j++)
01727 {
01728 work[j] = 1.;
01729
01730 for (octave_idx_type k = j; k >= 0; k--)
01731 {
01732 if (work[k] != 0.)
01733 {
01734 Complex tmp = work[k] / data(cidx(k+1)-1);
01735 work[k] = tmp;
01736 for (octave_idx_type i = cidx(k);
01737 i < cidx(k+1)-1; i++)
01738 {
01739 octave_idx_type iidx = ridx(i);
01740 work[iidx] = work[iidx] - tmp * data(i);
01741 }
01742 }
01743 }
01744 double atmp = 0;
01745 for (octave_idx_type i = 0; i < j+1; i++)
01746 {
01747 atmp += std::abs(work[i]);
01748 work[i] = 0.;
01749 }
01750 if (atmp > ainvnorm)
01751 ainvnorm = atmp;
01752 }
01753 rcond = 1. / ainvnorm / anorm;
01754 }
01755 }
01756
01757 triangular_error:
01758 if (err != 0)
01759 {
01760 if (sing_handler)
01761 {
01762 sing_handler (rcond);
01763 mattype.mark_as_rectangular ();
01764 }
01765 else
01766 (*current_liboctave_error_handler)
01767 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
01768 rcond);
01769 }
01770
01771 volatile double rcond_plus_one = rcond + 1.0;
01772
01773 if (rcond_plus_one == 1.0 || xisnan (rcond))
01774 {
01775 err = -2;
01776
01777 if (sing_handler)
01778 {
01779 sing_handler (rcond);
01780 mattype.mark_as_rectangular ();
01781 }
01782 else
01783 (*current_liboctave_error_handler)
01784 ("matrix singular to machine precision, rcond = %g",
01785 rcond);
01786 }
01787 }
01788 else
01789 (*current_liboctave_error_handler) ("incorrect matrix type");
01790 }
01791
01792 return retval;
01793 }
01794
01795 SparseComplexMatrix
01796 SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseMatrix& b,
01797 octave_idx_type& err, double& rcond,
01798 solve_singularity_handler sing_handler,
01799 bool calc_cond) const
01800 {
01801 SparseComplexMatrix retval;
01802
01803 octave_idx_type nr = rows ();
01804 octave_idx_type nc = cols ();
01805 octave_idx_type nm = (nc > nr ? nc : nr);
01806 err = 0;
01807
01808 if (nr != b.rows ())
01809 (*current_liboctave_error_handler)
01810 ("matrix dimension mismatch solution of linear equations");
01811 else if (nr == 0 || nc == 0 || b.cols () == 0)
01812 retval = SparseComplexMatrix (nc, b.cols ());
01813 else
01814 {
01815
01816 int typ = mattype.type ();
01817 mattype.info ();
01818
01819 if (typ == MatrixType::Permuted_Upper ||
01820 typ == MatrixType::Upper)
01821 {
01822 double anorm = 0.;
01823 double ainvnorm = 0.;
01824 rcond = 1.;
01825
01826 if (calc_cond)
01827 {
01828
01829 for (octave_idx_type j = 0; j < nc; j++)
01830 {
01831 double atmp = 0.;
01832 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
01833 atmp += std::abs(data(i));
01834 if (atmp > anorm)
01835 anorm = atmp;
01836 }
01837 }
01838
01839 octave_idx_type b_nc = b.cols ();
01840 octave_idx_type b_nz = b.nnz ();
01841 retval = SparseComplexMatrix (nc, b_nc, b_nz);
01842 retval.xcidx(0) = 0;
01843 octave_idx_type ii = 0;
01844 octave_idx_type x_nz = b_nz;
01845
01846 if (typ == MatrixType::Permuted_Upper)
01847 {
01848 octave_idx_type *perm = mattype.triangular_perm ();
01849 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
01850
01851 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
01852 for (octave_idx_type i = 0; i < nc; i++)
01853 rperm[perm[i]] = i;
01854
01855 for (octave_idx_type j = 0; j < b_nc; j++)
01856 {
01857 for (octave_idx_type i = 0; i < nm; i++)
01858 work[i] = 0.;
01859 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
01860 work[b.ridx(i)] = b.data(i);
01861
01862 for (octave_idx_type k = nc-1; k >= 0; k--)
01863 {
01864 octave_idx_type kidx = perm[k];
01865
01866 if (work[k] != 0.)
01867 {
01868 if (ridx(cidx(kidx+1)-1) != k ||
01869 data(cidx(kidx+1)-1) == 0.)
01870 {
01871 err = -2;
01872 goto triangular_error;
01873 }
01874
01875 Complex tmp = work[k] / data(cidx(kidx+1)-1);
01876 work[k] = tmp;
01877 for (octave_idx_type i = cidx(kidx);
01878 i < cidx(kidx+1)-1; i++)
01879 {
01880 octave_idx_type iidx = ridx(i);
01881 work[iidx] = work[iidx] - tmp * data(i);
01882 }
01883 }
01884 }
01885
01886
01887
01888 octave_idx_type new_nnz = 0;
01889 for (octave_idx_type i = 0; i < nc; i++)
01890 if (work[i] != 0.)
01891 new_nnz++;
01892
01893 if (ii + new_nnz > x_nz)
01894 {
01895
01896 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
01897 retval.change_capacity (sz);
01898 x_nz = sz;
01899 }
01900
01901 for (octave_idx_type i = 0; i < nc; i++)
01902 if (work[rperm[i]] != 0.)
01903 {
01904 retval.xridx(ii) = i;
01905 retval.xdata(ii++) = work[rperm[i]];
01906 }
01907 retval.xcidx(j+1) = ii;
01908 }
01909
01910 retval.maybe_compress ();
01911
01912 if (calc_cond)
01913 {
01914
01915 for (octave_idx_type i = 0; i < nm; i++)
01916 work[i] = 0.;
01917
01918 for (octave_idx_type j = 0; j < nr; j++)
01919 {
01920 work[j] = 1.;
01921
01922 for (octave_idx_type k = j; k >= 0; k--)
01923 {
01924 octave_idx_type iidx = perm[k];
01925
01926 if (work[k] != 0.)
01927 {
01928 Complex tmp = work[k] / data(cidx(iidx+1)-1);
01929 work[k] = tmp;
01930 for (octave_idx_type i = cidx(iidx);
01931 i < cidx(iidx+1)-1; i++)
01932 {
01933 octave_idx_type idx2 = ridx(i);
01934 work[idx2] = work[idx2] - tmp * data(i);
01935 }
01936 }
01937 }
01938 double atmp = 0;
01939 for (octave_idx_type i = 0; i < j+1; i++)
01940 {
01941 atmp += std::abs(work[i]);
01942 work[i] = 0.;
01943 }
01944 if (atmp > ainvnorm)
01945 ainvnorm = atmp;
01946 }
01947 rcond = 1. / ainvnorm / anorm;
01948 }
01949 }
01950 else
01951 {
01952 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
01953
01954 for (octave_idx_type j = 0; j < b_nc; j++)
01955 {
01956 for (octave_idx_type i = 0; i < nm; i++)
01957 work[i] = 0.;
01958 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
01959 work[b.ridx(i)] = b.data(i);
01960
01961 for (octave_idx_type k = nc-1; k >= 0; k--)
01962 {
01963 if (work[k] != 0.)
01964 {
01965 if (ridx(cidx(k+1)-1) != k ||
01966 data(cidx(k+1)-1) == 0.)
01967 {
01968 err = -2;
01969 goto triangular_error;
01970 }
01971
01972 Complex tmp = work[k] / data(cidx(k+1)-1);
01973 work[k] = tmp;
01974 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
01975 {
01976 octave_idx_type iidx = ridx(i);
01977 work[iidx] = work[iidx] - tmp * data(i);
01978 }
01979 }
01980 }
01981
01982
01983
01984 octave_idx_type new_nnz = 0;
01985 for (octave_idx_type i = 0; i < nc; i++)
01986 if (work[i] != 0.)
01987 new_nnz++;
01988
01989 if (ii + new_nnz > x_nz)
01990 {
01991
01992 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
01993 retval.change_capacity (sz);
01994 x_nz = sz;
01995 }
01996
01997 for (octave_idx_type i = 0; i < nc; i++)
01998 if (work[i] != 0.)
01999 {
02000 retval.xridx(ii) = i;
02001 retval.xdata(ii++) = work[i];
02002 }
02003 retval.xcidx(j+1) = ii;
02004 }
02005
02006 retval.maybe_compress ();
02007
02008 if (calc_cond)
02009 {
02010
02011 for (octave_idx_type i = 0; i < nm; i++)
02012 work[i] = 0.;
02013
02014 for (octave_idx_type j = 0; j < nr; j++)
02015 {
02016 work[j] = 1.;
02017
02018 for (octave_idx_type k = j; k >= 0; k--)
02019 {
02020 if (work[k] != 0.)
02021 {
02022 Complex tmp = work[k] / data(cidx(k+1)-1);
02023 work[k] = tmp;
02024 for (octave_idx_type i = cidx(k);
02025 i < cidx(k+1)-1; i++)
02026 {
02027 octave_idx_type iidx = ridx(i);
02028 work[iidx] = work[iidx] - tmp * data(i);
02029 }
02030 }
02031 }
02032 double atmp = 0;
02033 for (octave_idx_type i = 0; i < j+1; i++)
02034 {
02035 atmp += std::abs(work[i]);
02036 work[i] = 0.;
02037 }
02038 if (atmp > ainvnorm)
02039 ainvnorm = atmp;
02040 }
02041 rcond = 1. / ainvnorm / anorm;
02042 }
02043 }
02044
02045 triangular_error:
02046 if (err != 0)
02047 {
02048 if (sing_handler)
02049 {
02050 sing_handler (rcond);
02051 mattype.mark_as_rectangular ();
02052 }
02053 else
02054 (*current_liboctave_error_handler)
02055 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
02056 rcond);
02057 }
02058
02059 volatile double rcond_plus_one = rcond + 1.0;
02060
02061 if (rcond_plus_one == 1.0 || xisnan (rcond))
02062 {
02063 err = -2;
02064
02065 if (sing_handler)
02066 {
02067 sing_handler (rcond);
02068 mattype.mark_as_rectangular ();
02069 }
02070 else
02071 (*current_liboctave_error_handler)
02072 ("matrix singular to machine precision, rcond = %g",
02073 rcond);
02074 }
02075 }
02076 else
02077 (*current_liboctave_error_handler) ("incorrect matrix type");
02078 }
02079 return retval;
02080 }
02081
02082 ComplexMatrix
02083 SparseComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
02084 octave_idx_type& err, double& rcond,
02085 solve_singularity_handler sing_handler,
02086 bool calc_cond) const
02087 {
02088 ComplexMatrix retval;
02089
02090 octave_idx_type nr = rows ();
02091 octave_idx_type nc = cols ();
02092 octave_idx_type nm = (nc > nr ? nc : nr);
02093 err = 0;
02094
02095 if (nr != b.rows ())
02096 (*current_liboctave_error_handler)
02097 ("matrix dimension mismatch solution of linear equations");
02098 else if (nr == 0 || nc == 0 || b.cols () == 0)
02099 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
02100 else
02101 {
02102
02103 int typ = mattype.type ();
02104 mattype.info ();
02105
02106 if (typ == MatrixType::Permuted_Upper ||
02107 typ == MatrixType::Upper)
02108 {
02109 double anorm = 0.;
02110 double ainvnorm = 0.;
02111 octave_idx_type b_nc = b.cols ();
02112 rcond = 1.;
02113
02114 if (calc_cond)
02115 {
02116
02117 for (octave_idx_type j = 0; j < nc; j++)
02118 {
02119 double atmp = 0.;
02120 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
02121 atmp += std::abs(data(i));
02122 if (atmp > anorm)
02123 anorm = atmp;
02124 }
02125 }
02126
02127 if (typ == MatrixType::Permuted_Upper)
02128 {
02129 retval.resize (nc, b_nc);
02130 octave_idx_type *perm = mattype.triangular_perm ();
02131 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
02132
02133 for (octave_idx_type j = 0; j < b_nc; j++)
02134 {
02135 for (octave_idx_type i = 0; i < nr; i++)
02136 work[i] = b(i,j);
02137 for (octave_idx_type i = nr; i < nc; i++)
02138 work[i] = 0.;
02139
02140 for (octave_idx_type k = nc-1; k >= 0; k--)
02141 {
02142 octave_idx_type kidx = perm[k];
02143
02144 if (work[k] != 0.)
02145 {
02146 if (ridx(cidx(kidx+1)-1) != k ||
02147 data(cidx(kidx+1)-1) == 0.)
02148 {
02149 err = -2;
02150 goto triangular_error;
02151 }
02152
02153 Complex tmp = work[k] / data(cidx(kidx+1)-1);
02154 work[k] = tmp;
02155 for (octave_idx_type i = cidx(kidx);
02156 i < cidx(kidx+1)-1; i++)
02157 {
02158 octave_idx_type iidx = ridx(i);
02159 work[iidx] = work[iidx] - tmp * data(i);
02160 }
02161 }
02162 }
02163
02164 for (octave_idx_type i = 0; i < nc; i++)
02165 retval (perm[i], j) = work[i];
02166 }
02167
02168 if (calc_cond)
02169 {
02170
02171 for (octave_idx_type i = 0; i < nm; i++)
02172 work[i] = 0.;
02173
02174 for (octave_idx_type j = 0; j < nr; j++)
02175 {
02176 work[j] = 1.;
02177
02178 for (octave_idx_type k = j; k >= 0; k--)
02179 {
02180 octave_idx_type iidx = perm[k];
02181
02182 if (work[k] != 0.)
02183 {
02184 Complex tmp = work[k] / data(cidx(iidx+1)-1);
02185 work[k] = tmp;
02186 for (octave_idx_type i = cidx(iidx);
02187 i < cidx(iidx+1)-1; i++)
02188 {
02189 octave_idx_type idx2 = ridx(i);
02190 work[idx2] = work[idx2] - tmp * data(i);
02191 }
02192 }
02193 }
02194 double atmp = 0;
02195 for (octave_idx_type i = 0; i < j+1; i++)
02196 {
02197 atmp += std::abs(work[i]);
02198 work[i] = 0.;
02199 }
02200 if (atmp > ainvnorm)
02201 ainvnorm = atmp;
02202 }
02203 rcond = 1. / ainvnorm / anorm;
02204 }
02205 }
02206 else
02207 {
02208 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
02209 retval.resize (nc, b_nc);
02210
02211 for (octave_idx_type j = 0; j < b_nc; j++)
02212 {
02213 for (octave_idx_type i = 0; i < nr; i++)
02214 work[i] = b(i,j);
02215 for (octave_idx_type i = nr; i < nc; i++)
02216 work[i] = 0.;
02217
02218 for (octave_idx_type k = nc-1; k >= 0; k--)
02219 {
02220 if (work[k] != 0.)
02221 {
02222 if (ridx(cidx(k+1)-1) != k ||
02223 data(cidx(k+1)-1) == 0.)
02224 {
02225 err = -2;
02226 goto triangular_error;
02227 }
02228
02229 Complex tmp = work[k] / data(cidx(k+1)-1);
02230 work[k] = tmp;
02231 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
02232 {
02233 octave_idx_type iidx = ridx(i);
02234 work[iidx] = work[iidx] - tmp * data(i);
02235 }
02236 }
02237 }
02238
02239 for (octave_idx_type i = 0; i < nc; i++)
02240 retval.xelem (i, j) = work[i];
02241 }
02242
02243 if (calc_cond)
02244 {
02245
02246 for (octave_idx_type i = 0; i < nm; i++)
02247 work[i] = 0.;
02248
02249 for (octave_idx_type j = 0; j < nr; j++)
02250 {
02251 work[j] = 1.;
02252
02253 for (octave_idx_type k = j; k >= 0; k--)
02254 {
02255 if (work[k] != 0.)
02256 {
02257 Complex tmp = work[k] / data(cidx(k+1)-1);
02258 work[k] = tmp;
02259 for (octave_idx_type i = cidx(k);
02260 i < cidx(k+1)-1; i++)
02261 {
02262 octave_idx_type iidx = ridx(i);
02263 work[iidx] = work[iidx] - tmp * data(i);
02264 }
02265 }
02266 }
02267 double atmp = 0;
02268 for (octave_idx_type i = 0; i < j+1; i++)
02269 {
02270 atmp += std::abs(work[i]);
02271 work[i] = 0.;
02272 }
02273 if (atmp > ainvnorm)
02274 ainvnorm = atmp;
02275 }
02276 rcond = 1. / ainvnorm / anorm;
02277 }
02278 }
02279
02280 triangular_error:
02281 if (err != 0)
02282 {
02283 if (sing_handler)
02284 {
02285 sing_handler (rcond);
02286 mattype.mark_as_rectangular ();
02287 }
02288 else
02289 (*current_liboctave_error_handler)
02290 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
02291 rcond);
02292 }
02293
02294 volatile double rcond_plus_one = rcond + 1.0;
02295
02296 if (rcond_plus_one == 1.0 || xisnan (rcond))
02297 {
02298 err = -2;
02299
02300 if (sing_handler)
02301 {
02302 sing_handler (rcond);
02303 mattype.mark_as_rectangular ();
02304 }
02305 else
02306 (*current_liboctave_error_handler)
02307 ("matrix singular to machine precision, rcond = %g",
02308 rcond);
02309 }
02310 }
02311 else
02312 (*current_liboctave_error_handler) ("incorrect matrix type");
02313 }
02314
02315 return retval;
02316 }
02317
02318 SparseComplexMatrix
02319 SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseComplexMatrix& b,
02320 octave_idx_type& err, double& rcond,
02321 solve_singularity_handler sing_handler,
02322 bool calc_cond) const
02323 {
02324 SparseComplexMatrix retval;
02325
02326 octave_idx_type nr = rows ();
02327 octave_idx_type nc = cols ();
02328 octave_idx_type nm = (nc > nr ? nc : nr);
02329 err = 0;
02330
02331 if (nr != b.rows ())
02332 (*current_liboctave_error_handler)
02333 ("matrix dimension mismatch solution of linear equations");
02334 else if (nr == 0 || nc == 0 || b.cols () == 0)
02335 retval = SparseComplexMatrix (nc, b.cols ());
02336 else
02337 {
02338
02339 int typ = mattype.type ();
02340 mattype.info ();
02341
02342 if (typ == MatrixType::Permuted_Upper ||
02343 typ == MatrixType::Upper)
02344 {
02345 double anorm = 0.;
02346 double ainvnorm = 0.;
02347 rcond = 1.;
02348
02349 if (calc_cond)
02350 {
02351
02352 for (octave_idx_type j = 0; j < nc; j++)
02353 {
02354 double atmp = 0.;
02355 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
02356 atmp += std::abs(data(i));
02357 if (atmp > anorm)
02358 anorm = atmp;
02359 }
02360 }
02361
02362 octave_idx_type b_nc = b.cols ();
02363 octave_idx_type b_nz = b.nnz ();
02364 retval = SparseComplexMatrix (nc, b_nc, b_nz);
02365 retval.xcidx(0) = 0;
02366 octave_idx_type ii = 0;
02367 octave_idx_type x_nz = b_nz;
02368
02369 if (typ == MatrixType::Permuted_Upper)
02370 {
02371 octave_idx_type *perm = mattype.triangular_perm ();
02372 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
02373
02374 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
02375 for (octave_idx_type i = 0; i < nc; i++)
02376 rperm[perm[i]] = i;
02377
02378 for (octave_idx_type j = 0; j < b_nc; j++)
02379 {
02380 for (octave_idx_type i = 0; i < nm; i++)
02381 work[i] = 0.;
02382 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
02383 work[b.ridx(i)] = b.data(i);
02384
02385 for (octave_idx_type k = nc-1; k >= 0; k--)
02386 {
02387 octave_idx_type kidx = perm[k];
02388
02389 if (work[k] != 0.)
02390 {
02391 if (ridx(cidx(kidx+1)-1) != k ||
02392 data(cidx(kidx+1)-1) == 0.)
02393 {
02394 err = -2;
02395 goto triangular_error;
02396 }
02397
02398 Complex tmp = work[k] / data(cidx(kidx+1)-1);
02399 work[k] = tmp;
02400 for (octave_idx_type i = cidx(kidx);
02401 i < cidx(kidx+1)-1; i++)
02402 {
02403 octave_idx_type iidx = ridx(i);
02404 work[iidx] = work[iidx] - tmp * data(i);
02405 }
02406 }
02407 }
02408
02409
02410
02411 octave_idx_type new_nnz = 0;
02412 for (octave_idx_type i = 0; i < nc; i++)
02413 if (work[i] != 0.)
02414 new_nnz++;
02415
02416 if (ii + new_nnz > x_nz)
02417 {
02418
02419 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
02420 retval.change_capacity (sz);
02421 x_nz = sz;
02422 }
02423
02424 for (octave_idx_type i = 0; i < nc; i++)
02425 if (work[rperm[i]] != 0.)
02426 {
02427 retval.xridx(ii) = i;
02428 retval.xdata(ii++) = work[rperm[i]];
02429 }
02430 retval.xcidx(j+1) = ii;
02431 }
02432
02433 retval.maybe_compress ();
02434
02435 if (calc_cond)
02436 {
02437
02438 for (octave_idx_type i = 0; i < nm; i++)
02439 work[i] = 0.;
02440
02441 for (octave_idx_type j = 0; j < nr; j++)
02442 {
02443 work[j] = 1.;
02444
02445 for (octave_idx_type k = j; k >= 0; k--)
02446 {
02447 octave_idx_type iidx = perm[k];
02448
02449 if (work[k] != 0.)
02450 {
02451 Complex tmp = work[k] / data(cidx(iidx+1)-1);
02452 work[k] = tmp;
02453 for (octave_idx_type i = cidx(iidx);
02454 i < cidx(iidx+1)-1; i++)
02455 {
02456 octave_idx_type idx2 = ridx(i);
02457 work[idx2] = work[idx2] - tmp * data(i);
02458 }
02459 }
02460 }
02461 double atmp = 0;
02462 for (octave_idx_type i = 0; i < j+1; i++)
02463 {
02464 atmp += std::abs(work[i]);
02465 work[i] = 0.;
02466 }
02467 if (atmp > ainvnorm)
02468 ainvnorm = atmp;
02469 }
02470 rcond = 1. / ainvnorm / anorm;
02471 }
02472 }
02473 else
02474 {
02475 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
02476
02477 for (octave_idx_type j = 0; j < b_nc; j++)
02478 {
02479 for (octave_idx_type i = 0; i < nm; i++)
02480 work[i] = 0.;
02481 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
02482 work[b.ridx(i)] = b.data(i);
02483
02484 for (octave_idx_type k = nr-1; k >= 0; k--)
02485 {
02486 if (work[k] != 0.)
02487 {
02488 if (ridx(cidx(k+1)-1) != k ||
02489 data(cidx(k+1)-1) == 0.)
02490 {
02491 err = -2;
02492 goto triangular_error;
02493 }
02494
02495 Complex tmp = work[k] / data(cidx(k+1)-1);
02496 work[k] = tmp;
02497 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
02498 {
02499 octave_idx_type iidx = ridx(i);
02500 work[iidx] = work[iidx] - tmp * data(i);
02501 }
02502 }
02503 }
02504
02505
02506
02507 octave_idx_type new_nnz = 0;
02508 for (octave_idx_type i = 0; i < nc; i++)
02509 if (work[i] != 0.)
02510 new_nnz++;
02511
02512 if (ii + new_nnz > x_nz)
02513 {
02514
02515 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
02516 retval.change_capacity (sz);
02517 x_nz = sz;
02518 }
02519
02520 for (octave_idx_type i = 0; i < nc; i++)
02521 if (work[i] != 0.)
02522 {
02523 retval.xridx(ii) = i;
02524 retval.xdata(ii++) = work[i];
02525 }
02526 retval.xcidx(j+1) = ii;
02527 }
02528
02529 retval.maybe_compress ();
02530
02531 if (calc_cond)
02532 {
02533
02534 for (octave_idx_type i = 0; i < nm; i++)
02535 work[i] = 0.;
02536
02537 for (octave_idx_type j = 0; j < nr; j++)
02538 {
02539 work[j] = 1.;
02540
02541 for (octave_idx_type k = j; k >= 0; k--)
02542 {
02543 if (work[k] != 0.)
02544 {
02545 Complex tmp = work[k] / data(cidx(k+1)-1);
02546 work[k] = tmp;
02547 for (octave_idx_type i = cidx(k);
02548 i < cidx(k+1)-1; i++)
02549 {
02550 octave_idx_type iidx = ridx(i);
02551 work[iidx] = work[iidx] - tmp * data(i);
02552 }
02553 }
02554 }
02555 double atmp = 0;
02556 for (octave_idx_type i = 0; i < j+1; i++)
02557 {
02558 atmp += std::abs(work[i]);
02559 work[i] = 0.;
02560 }
02561 if (atmp > ainvnorm)
02562 ainvnorm = atmp;
02563 }
02564 rcond = 1. / ainvnorm / anorm;
02565 }
02566 }
02567
02568 triangular_error:
02569 if (err != 0)
02570 {
02571 if (sing_handler)
02572 {
02573 sing_handler (rcond);
02574 mattype.mark_as_rectangular ();
02575 }
02576 else
02577 (*current_liboctave_error_handler)
02578 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
02579 rcond);
02580 }
02581
02582 volatile double rcond_plus_one = rcond + 1.0;
02583
02584 if (rcond_plus_one == 1.0 || xisnan (rcond))
02585 {
02586 err = -2;
02587
02588 if (sing_handler)
02589 {
02590 sing_handler (rcond);
02591 mattype.mark_as_rectangular ();
02592 }
02593 else
02594 (*current_liboctave_error_handler)
02595 ("matrix singular to machine precision, rcond = %g",
02596 rcond);
02597 }
02598 }
02599 else
02600 (*current_liboctave_error_handler) ("incorrect matrix type");
02601 }
02602
02603 return retval;
02604 }
02605
02606 ComplexMatrix
02607 SparseComplexMatrix::ltsolve (MatrixType &mattype, const Matrix& b,
02608 octave_idx_type& err, double& rcond,
02609 solve_singularity_handler sing_handler,
02610 bool calc_cond) const
02611 {
02612 ComplexMatrix retval;
02613
02614 octave_idx_type nr = rows ();
02615 octave_idx_type nc = cols ();
02616 octave_idx_type nm = (nc > nr ? nc : nr);
02617 err = 0;
02618
02619 if (nr != b.rows ())
02620 (*current_liboctave_error_handler)
02621 ("matrix dimension mismatch solution of linear equations");
02622 else if (nr == 0 || nc == 0 || b.cols () == 0)
02623 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
02624 else
02625 {
02626
02627 int typ = mattype.type ();
02628 mattype.info ();
02629
02630 if (typ == MatrixType::Permuted_Lower ||
02631 typ == MatrixType::Lower)
02632 {
02633 double anorm = 0.;
02634 double ainvnorm = 0.;
02635 octave_idx_type b_nc = b.cols ();
02636 rcond = 1.;
02637
02638 if (calc_cond)
02639 {
02640
02641 for (octave_idx_type j = 0; j < nc; j++)
02642 {
02643 double atmp = 0.;
02644 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
02645 atmp += std::abs(data(i));
02646 if (atmp > anorm)
02647 anorm = atmp;
02648 }
02649 }
02650
02651 if (typ == MatrixType::Permuted_Lower)
02652 {
02653 retval.resize (nc, b_nc);
02654 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
02655 octave_idx_type *perm = mattype.triangular_perm ();
02656
02657 for (octave_idx_type j = 0; j < b_nc; j++)
02658 {
02659 for (octave_idx_type i = 0; i < nm; i++)
02660 work[i] = 0.;
02661 for (octave_idx_type i = 0; i < nr; i++)
02662 work[perm[i]] = b(i,j);
02663
02664 for (octave_idx_type k = 0; k < nc; k++)
02665 {
02666 if (work[k] != 0.)
02667 {
02668 octave_idx_type minr = nr;
02669 octave_idx_type mini = 0;
02670
02671 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
02672 if (perm[ridx(i)] < minr)
02673 {
02674 minr = perm[ridx(i)];
02675 mini = i;
02676 }
02677
02678 if (minr != k || data (mini) == 0.)
02679 {
02680 err = -2;
02681 goto triangular_error;
02682 }
02683
02684 Complex tmp = work[k] / data(mini);
02685 work[k] = tmp;
02686 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
02687 {
02688 if (i == mini)
02689 continue;
02690
02691 octave_idx_type iidx = perm[ridx(i)];
02692 work[iidx] = work[iidx] - tmp * data(i);
02693 }
02694 }
02695 }
02696
02697 for (octave_idx_type i = 0; i < nc; i++)
02698 retval (i, j) = work[i];
02699 }
02700
02701 if (calc_cond)
02702 {
02703
02704 for (octave_idx_type i = 0; i < nm; i++)
02705 work[i] = 0.;
02706
02707 for (octave_idx_type j = 0; j < nr; j++)
02708 {
02709 work[j] = 1.;
02710
02711 for (octave_idx_type k = 0; k < nc; k++)
02712 {
02713 if (work[k] != 0.)
02714 {
02715 octave_idx_type minr = nr;
02716 octave_idx_type mini = 0;
02717
02718 for (octave_idx_type i = cidx(k);
02719 i < cidx(k+1); i++)
02720 if (perm[ridx(i)] < minr)
02721 {
02722 minr = perm[ridx(i)];
02723 mini = i;
02724 }
02725
02726 Complex tmp = work[k] / data(mini);
02727 work[k] = tmp;
02728 for (octave_idx_type i = cidx(k);
02729 i < cidx(k+1); i++)
02730 {
02731 if (i == mini)
02732 continue;
02733
02734 octave_idx_type iidx = perm[ridx(i)];
02735 work[iidx] = work[iidx] - tmp * data(i);
02736 }
02737 }
02738 }
02739
02740 double atmp = 0;
02741 for (octave_idx_type i = j; i < nc; i++)
02742 {
02743 atmp += std::abs(work[i]);
02744 work[i] = 0.;
02745 }
02746 if (atmp > ainvnorm)
02747 ainvnorm = atmp;
02748 }
02749 rcond = 1. / ainvnorm / anorm;
02750 }
02751 }
02752 else
02753 {
02754 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
02755 retval.resize (nc, b_nc, 0.);
02756
02757 for (octave_idx_type j = 0; j < b_nc; j++)
02758 {
02759 for (octave_idx_type i = 0; i < nr; i++)
02760 work[i] = b(i,j);
02761 for (octave_idx_type i = nr; i < nc; i++)
02762 work[i] = 0.;
02763 for (octave_idx_type k = 0; k < nc; k++)
02764 {
02765 if (work[k] != 0.)
02766 {
02767 if (ridx(cidx(k)) != k ||
02768 data(cidx(k)) == 0.)
02769 {
02770 err = -2;
02771 goto triangular_error;
02772 }
02773
02774 Complex tmp = work[k] / data(cidx(k));
02775 work[k] = tmp;
02776 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
02777 {
02778 octave_idx_type iidx = ridx(i);
02779 work[iidx] = work[iidx] - tmp * data(i);
02780 }
02781 }
02782 }
02783 for (octave_idx_type i = 0; i < nc; i++)
02784 retval.xelem (i, j) = work[i];
02785 }
02786
02787 if (calc_cond)
02788 {
02789
02790 for (octave_idx_type i = 0; i < nm; i++)
02791 work[i] = 0.;
02792
02793 for (octave_idx_type j = 0; j < nr; j++)
02794 {
02795 work[j] = 1.;
02796
02797 for (octave_idx_type k = j; k < nc; k++)
02798 {
02799
02800 if (work[k] != 0.)
02801 {
02802 Complex tmp = work[k] / data(cidx(k));
02803 work[k] = tmp;
02804 for (octave_idx_type i = cidx(k)+1;
02805 i < cidx(k+1); i++)
02806 {
02807 octave_idx_type iidx = ridx(i);
02808 work[iidx] = work[iidx] - tmp * data(i);
02809 }
02810 }
02811 }
02812 double atmp = 0;
02813 for (octave_idx_type i = j; i < nc; i++)
02814 {
02815 atmp += std::abs(work[i]);
02816 work[i] = 0.;
02817 }
02818 if (atmp > ainvnorm)
02819 ainvnorm = atmp;
02820 }
02821 rcond = 1. / ainvnorm / anorm;
02822 }
02823 }
02824 triangular_error:
02825 if (err != 0)
02826 {
02827 if (sing_handler)
02828 {
02829 sing_handler (rcond);
02830 mattype.mark_as_rectangular ();
02831 }
02832 else
02833 (*current_liboctave_error_handler)
02834 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
02835 rcond);
02836 }
02837
02838 volatile double rcond_plus_one = rcond + 1.0;
02839
02840 if (rcond_plus_one == 1.0 || xisnan (rcond))
02841 {
02842 err = -2;
02843
02844 if (sing_handler)
02845 {
02846 sing_handler (rcond);
02847 mattype.mark_as_rectangular ();
02848 }
02849 else
02850 (*current_liboctave_error_handler)
02851 ("matrix singular to machine precision, rcond = %g",
02852 rcond);
02853 }
02854 }
02855 else
02856 (*current_liboctave_error_handler) ("incorrect matrix type");
02857 }
02858
02859 return retval;
02860 }
02861
02862 SparseComplexMatrix
02863 SparseComplexMatrix::ltsolve (MatrixType &mattype, const SparseMatrix& b,
02864 octave_idx_type& err, double& rcond,
02865 solve_singularity_handler sing_handler,
02866 bool calc_cond) const
02867 {
02868 SparseComplexMatrix retval;
02869
02870 octave_idx_type nr = rows ();
02871 octave_idx_type nc = cols ();
02872 octave_idx_type nm = (nc > nr ? nc : nr);
02873
02874 err = 0;
02875
02876 if (nr != b.rows ())
02877 (*current_liboctave_error_handler)
02878 ("matrix dimension mismatch solution of linear equations");
02879 else if (nr == 0 || nc == 0 || b.cols () == 0)
02880 retval = SparseComplexMatrix (nc, b.cols ());
02881 else
02882 {
02883
02884 int typ = mattype.type ();
02885 mattype.info ();
02886
02887 if (typ == MatrixType::Permuted_Lower ||
02888 typ == MatrixType::Lower)
02889 {
02890 double anorm = 0.;
02891 double ainvnorm = 0.;
02892 rcond = 1.;
02893
02894 if (calc_cond)
02895 {
02896
02897 for (octave_idx_type j = 0; j < nc; j++)
02898 {
02899 double atmp = 0.;
02900 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
02901 atmp += std::abs(data(i));
02902 if (atmp > anorm)
02903 anorm = atmp;
02904 }
02905 }
02906
02907 octave_idx_type b_nc = b.cols ();
02908 octave_idx_type b_nz = b.nnz ();
02909 retval = SparseComplexMatrix (nc, b_nc, b_nz);
02910 retval.xcidx(0) = 0;
02911 octave_idx_type ii = 0;
02912 octave_idx_type x_nz = b_nz;
02913
02914 if (typ == MatrixType::Permuted_Lower)
02915 {
02916 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
02917 octave_idx_type *perm = mattype.triangular_perm ();
02918
02919 for (octave_idx_type j = 0; j < b_nc; j++)
02920 {
02921 for (octave_idx_type i = 0; i < nm; i++)
02922 work[i] = 0.;
02923 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
02924 work[perm[b.ridx(i)]] = b.data(i);
02925
02926 for (octave_idx_type k = 0; k < nc; k++)
02927 {
02928 if (work[k] != 0.)
02929 {
02930 octave_idx_type minr = nr;
02931 octave_idx_type mini = 0;
02932
02933 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
02934 if (perm[ridx(i)] < minr)
02935 {
02936 minr = perm[ridx(i)];
02937 mini = i;
02938 }
02939
02940 if (minr != k || data (mini) == 0.)
02941 {
02942 err = -2;
02943 goto triangular_error;
02944 }
02945
02946 Complex tmp = work[k] / data(mini);
02947 work[k] = tmp;
02948 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
02949 {
02950 if (i == mini)
02951 continue;
02952
02953 octave_idx_type iidx = perm[ridx(i)];
02954 work[iidx] = work[iidx] - tmp * data(i);
02955 }
02956 }
02957 }
02958
02959
02960
02961 octave_idx_type new_nnz = 0;
02962 for (octave_idx_type i = 0; i < nc; i++)
02963 if (work[i] != 0.)
02964 new_nnz++;
02965
02966 if (ii + new_nnz > x_nz)
02967 {
02968
02969 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
02970 retval.change_capacity (sz);
02971 x_nz = sz;
02972 }
02973
02974 for (octave_idx_type i = 0; i < nc; i++)
02975 if (work[i] != 0.)
02976 {
02977 retval.xridx(ii) = i;
02978 retval.xdata(ii++) = work[i];
02979 }
02980 retval.xcidx(j+1) = ii;
02981 }
02982
02983 retval.maybe_compress ();
02984
02985 if (calc_cond)
02986 {
02987
02988 for (octave_idx_type i = 0; i < nm; i++)
02989 work[i] = 0.;
02990
02991 for (octave_idx_type j = 0; j < nr; j++)
02992 {
02993 work[j] = 1.;
02994
02995 for (octave_idx_type k = 0; k < nc; k++)
02996 {
02997 if (work[k] != 0.)
02998 {
02999 octave_idx_type minr = nr;
03000 octave_idx_type mini = 0;
03001
03002 for (octave_idx_type i = cidx(k);
03003 i < cidx(k+1); i++)
03004 if (perm[ridx(i)] < minr)
03005 {
03006 minr = perm[ridx(i)];
03007 mini = i;
03008 }
03009
03010 Complex tmp = work[k] / data(mini);
03011 work[k] = tmp;
03012 for (octave_idx_type i = cidx(k);
03013 i < cidx(k+1); i++)
03014 {
03015 if (i == mini)
03016 continue;
03017
03018 octave_idx_type iidx = perm[ridx(i)];
03019 work[iidx] = work[iidx] - tmp * data(i);
03020 }
03021 }
03022 }
03023
03024 double atmp = 0;
03025 for (octave_idx_type i = j; i < nc; i++)
03026 {
03027 atmp += std::abs(work[i]);
03028 work[i] = 0.;
03029 }
03030 if (atmp > ainvnorm)
03031 ainvnorm = atmp;
03032 }
03033 rcond = 1. / ainvnorm / anorm;
03034 }
03035 }
03036 else
03037 {
03038 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
03039
03040 for (octave_idx_type j = 0; j < b_nc; j++)
03041 {
03042 for (octave_idx_type i = 0; i < nm; i++)
03043 work[i] = 0.;
03044 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
03045 work[b.ridx(i)] = b.data(i);
03046
03047 for (octave_idx_type k = 0; k < nc; k++)
03048 {
03049 if (work[k] != 0.)
03050 {
03051 if (ridx(cidx(k)) != k ||
03052 data(cidx(k)) == 0.)
03053 {
03054 err = -2;
03055 goto triangular_error;
03056 }
03057
03058 Complex tmp = work[k] / data(cidx(k));
03059 work[k] = tmp;
03060 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
03061 {
03062 octave_idx_type iidx = ridx(i);
03063 work[iidx] = work[iidx] - tmp * data(i);
03064 }
03065 }
03066 }
03067
03068
03069
03070 octave_idx_type new_nnz = 0;
03071 for (octave_idx_type i = 0; i < nc; i++)
03072 if (work[i] != 0.)
03073 new_nnz++;
03074
03075 if (ii + new_nnz > x_nz)
03076 {
03077
03078 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
03079 retval.change_capacity (sz);
03080 x_nz = sz;
03081 }
03082
03083 for (octave_idx_type i = 0; i < nc; i++)
03084 if (work[i] != 0.)
03085 {
03086 retval.xridx(ii) = i;
03087 retval.xdata(ii++) = work[i];
03088 }
03089 retval.xcidx(j+1) = ii;
03090 }
03091
03092 retval.maybe_compress ();
03093
03094 if (calc_cond)
03095 {
03096
03097 for (octave_idx_type i = 0; i < nm; i++)
03098 work[i] = 0.;
03099
03100 for (octave_idx_type j = 0; j < nr; j++)
03101 {
03102 work[j] = 1.;
03103
03104 for (octave_idx_type k = j; k < nc; k++)
03105 {
03106
03107 if (work[k] != 0.)
03108 {
03109 Complex tmp = work[k] / data(cidx(k));
03110 work[k] = tmp;
03111 for (octave_idx_type i = cidx(k)+1;
03112 i < cidx(k+1); i++)
03113 {
03114 octave_idx_type iidx = ridx(i);
03115 work[iidx] = work[iidx] - tmp * data(i);
03116 }
03117 }
03118 }
03119 double atmp = 0;
03120 for (octave_idx_type i = j; i < nc; i++)
03121 {
03122 atmp += std::abs(work[i]);
03123 work[i] = 0.;
03124 }
03125 if (atmp > ainvnorm)
03126 ainvnorm = atmp;
03127 }
03128 rcond = 1. / ainvnorm / anorm;
03129 }
03130 }
03131
03132 triangular_error:
03133 if (err != 0)
03134 {
03135 if (sing_handler)
03136 {
03137 sing_handler (rcond);
03138 mattype.mark_as_rectangular ();
03139 }
03140 else
03141 (*current_liboctave_error_handler)
03142 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
03143 rcond);
03144 }
03145
03146 volatile double rcond_plus_one = rcond + 1.0;
03147
03148 if (rcond_plus_one == 1.0 || xisnan (rcond))
03149 {
03150 err = -2;
03151
03152 if (sing_handler)
03153 {
03154 sing_handler (rcond);
03155 mattype.mark_as_rectangular ();
03156 }
03157 else
03158 (*current_liboctave_error_handler)
03159 ("matrix singular to machine precision, rcond = %g",
03160 rcond);
03161 }
03162 }
03163 else
03164 (*current_liboctave_error_handler) ("incorrect matrix type");
03165 }
03166
03167 return retval;
03168 }
03169
03170 ComplexMatrix
03171 SparseComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
03172 octave_idx_type& err, double& rcond,
03173 solve_singularity_handler sing_handler,
03174 bool calc_cond) const
03175 {
03176 ComplexMatrix retval;
03177
03178 octave_idx_type nr = rows ();
03179 octave_idx_type nc = cols ();
03180 octave_idx_type nm = (nc > nr ? nc : nr);
03181 err = 0;
03182
03183 if (nr != b.rows ())
03184 (*current_liboctave_error_handler)
03185 ("matrix dimension mismatch solution of linear equations");
03186 else if (nr == 0 || nc == 0 || b.cols () == 0)
03187 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
03188 else
03189 {
03190
03191 int typ = mattype.type ();
03192 mattype.info ();
03193
03194 if (typ == MatrixType::Permuted_Lower ||
03195 typ == MatrixType::Lower)
03196 {
03197 double anorm = 0.;
03198 double ainvnorm = 0.;
03199 octave_idx_type b_nc = b.cols ();
03200 rcond = 1.;
03201
03202 if (calc_cond)
03203 {
03204
03205 for (octave_idx_type j = 0; j < nc; j++)
03206 {
03207 double atmp = 0.;
03208 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
03209 atmp += std::abs(data(i));
03210 if (atmp > anorm)
03211 anorm = atmp;
03212 }
03213 }
03214
03215 if (typ == MatrixType::Permuted_Lower)
03216 {
03217 retval.resize (nc, b_nc);
03218 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
03219 octave_idx_type *perm = mattype.triangular_perm ();
03220
03221 for (octave_idx_type j = 0; j < b_nc; j++)
03222 {
03223 for (octave_idx_type i = 0; i < nm; i++)
03224 work[i] = 0.;
03225 for (octave_idx_type i = 0; i < nr; i++)
03226 work[perm[i]] = b(i,j);
03227
03228 for (octave_idx_type k = 0; k < nc; k++)
03229 {
03230 if (work[k] != 0.)
03231 {
03232 octave_idx_type minr = nr;
03233 octave_idx_type mini = 0;
03234
03235 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
03236 if (perm[ridx(i)] < minr)
03237 {
03238 minr = perm[ridx(i)];
03239 mini = i;
03240 }
03241
03242 if (minr != k || data (mini) == 0.)
03243 {
03244 err = -2;
03245 goto triangular_error;
03246 }
03247
03248 Complex tmp = work[k] / data(mini);
03249 work[k] = tmp;
03250 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
03251 {
03252 if (i == mini)
03253 continue;
03254
03255 octave_idx_type iidx = perm[ridx(i)];
03256 work[iidx] = work[iidx] - tmp * data(i);
03257 }
03258 }
03259 }
03260
03261 for (octave_idx_type i = 0; i < nc; i++)
03262 retval (i, j) = work[i];
03263 }
03264
03265 if (calc_cond)
03266 {
03267
03268 for (octave_idx_type i = 0; i < nm; i++)
03269 work[i] = 0.;
03270
03271 for (octave_idx_type j = 0; j < nr; j++)
03272 {
03273 work[j] = 1.;
03274
03275 for (octave_idx_type k = 0; k < nc; k++)
03276 {
03277 if (work[k] != 0.)
03278 {
03279 octave_idx_type minr = nr;
03280 octave_idx_type mini = 0;
03281
03282 for (octave_idx_type i = cidx(k);
03283 i < cidx(k+1); i++)
03284 if (perm[ridx(i)] < minr)
03285 {
03286 minr = perm[ridx(i)];
03287 mini = i;
03288 }
03289
03290 Complex tmp = work[k] / data(mini);
03291 work[k] = tmp;
03292 for (octave_idx_type i = cidx(k);
03293 i < cidx(k+1); i++)
03294 {
03295 if (i == mini)
03296 continue;
03297
03298 octave_idx_type iidx = perm[ridx(i)];
03299 work[iidx] = work[iidx] - tmp * data(i);
03300 }
03301 }
03302 }
03303
03304 double atmp = 0;
03305 for (octave_idx_type i = j; i < nc; i++)
03306 {
03307 atmp += std::abs(work[i]);
03308 work[i] = 0.;
03309 }
03310 if (atmp > ainvnorm)
03311 ainvnorm = atmp;
03312 }
03313 rcond = 1. / ainvnorm / anorm;
03314 }
03315 }
03316 else
03317 {
03318 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
03319 retval.resize (nc, b_nc, 0.);
03320
03321
03322 for (octave_idx_type j = 0; j < b_nc; j++)
03323 {
03324 for (octave_idx_type i = 0; i < nr; i++)
03325 work[i] = b(i,j);
03326 for (octave_idx_type i = nr; i < nc; i++)
03327 work[i] = 0.;
03328
03329 for (octave_idx_type k = 0; k < nc; k++)
03330 {
03331 if (work[k] != 0.)
03332 {
03333 if (ridx(cidx(k)) != k ||
03334 data(cidx(k)) == 0.)
03335 {
03336 err = -2;
03337 goto triangular_error;
03338 }
03339
03340 Complex tmp = work[k] / data(cidx(k));
03341 work[k] = tmp;
03342 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
03343 {
03344 octave_idx_type iidx = ridx(i);
03345 work[iidx] = work[iidx] - tmp * data(i);
03346 }
03347 }
03348 }
03349
03350 for (octave_idx_type i = 0; i < nc; i++)
03351 retval.xelem (i, j) = work[i];
03352 }
03353
03354 if (calc_cond)
03355 {
03356
03357 for (octave_idx_type i = 0; i < nm; i++)
03358 work[i] = 0.;
03359
03360 for (octave_idx_type j = 0; j < nr; j++)
03361 {
03362 work[j] = 1.;
03363
03364 for (octave_idx_type k = j; k < nc; k++)
03365 {
03366
03367 if (work[k] != 0.)
03368 {
03369 Complex tmp = work[k] / data(cidx(k));
03370 work[k] = tmp;
03371 for (octave_idx_type i = cidx(k)+1;
03372 i < cidx(k+1); i++)
03373 {
03374 octave_idx_type iidx = ridx(i);
03375 work[iidx] = work[iidx] - tmp * data(i);
03376 }
03377 }
03378 }
03379 double atmp = 0;
03380 for (octave_idx_type i = j; i < nc; i++)
03381 {
03382 atmp += std::abs(work[i]);
03383 work[i] = 0.;
03384 }
03385 if (atmp > ainvnorm)
03386 ainvnorm = atmp;
03387 }
03388 rcond = 1. / ainvnorm / anorm;
03389 }
03390 }
03391
03392 triangular_error:
03393 if (err != 0)
03394 {
03395 if (sing_handler)
03396 {
03397 sing_handler (rcond);
03398 mattype.mark_as_rectangular ();
03399 }
03400 else
03401 (*current_liboctave_error_handler)
03402 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
03403 rcond);
03404 }
03405
03406 volatile double rcond_plus_one = rcond + 1.0;
03407
03408 if (rcond_plus_one == 1.0 || xisnan (rcond))
03409 {
03410 err = -2;
03411
03412 if (sing_handler)
03413 {
03414 sing_handler (rcond);
03415 mattype.mark_as_rectangular ();
03416 }
03417 else
03418 (*current_liboctave_error_handler)
03419 ("matrix singular to machine precision, rcond = %g",
03420 rcond);
03421 }
03422 }
03423 else
03424 (*current_liboctave_error_handler) ("incorrect matrix type");
03425 }
03426
03427 return retval;
03428 }
03429
03430 SparseComplexMatrix
03431 SparseComplexMatrix::ltsolve (MatrixType &mattype, const SparseComplexMatrix& b,
03432 octave_idx_type& err, double& rcond,
03433 solve_singularity_handler sing_handler,
03434 bool calc_cond) const
03435 {
03436 SparseComplexMatrix retval;
03437
03438 octave_idx_type nr = rows ();
03439 octave_idx_type nc = cols ();
03440 octave_idx_type nm = (nc > nr ? nc : nr);
03441 err = 0;
03442
03443 if (nr != b.rows ())
03444 (*current_liboctave_error_handler)
03445 ("matrix dimension mismatch solution of linear equations");
03446 else if (nr == 0 || nc == 0 || b.cols () == 0)
03447 retval = SparseComplexMatrix (nc, b.cols ());
03448 else
03449 {
03450
03451 int typ = mattype.type ();
03452 mattype.info ();
03453
03454 if (typ == MatrixType::Permuted_Lower ||
03455 typ == MatrixType::Lower)
03456 {
03457 double anorm = 0.;
03458 double ainvnorm = 0.;
03459 rcond = 1.;
03460
03461 if (calc_cond)
03462 {
03463
03464 for (octave_idx_type j = 0; j < nc; j++)
03465 {
03466 double atmp = 0.;
03467 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
03468 atmp += std::abs(data(i));
03469 if (atmp > anorm)
03470 anorm = atmp;
03471 }
03472 }
03473
03474 octave_idx_type b_nc = b.cols ();
03475 octave_idx_type b_nz = b.nnz ();
03476 retval = SparseComplexMatrix (nc, b_nc, b_nz);
03477 retval.xcidx(0) = 0;
03478 octave_idx_type ii = 0;
03479 octave_idx_type x_nz = b_nz;
03480
03481 if (typ == MatrixType::Permuted_Lower)
03482 {
03483 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
03484 octave_idx_type *perm = mattype.triangular_perm ();
03485
03486 for (octave_idx_type j = 0; j < b_nc; j++)
03487 {
03488 for (octave_idx_type i = 0; i < nm; i++)
03489 work[i] = 0.;
03490 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
03491 work[perm[b.ridx(i)]] = b.data(i);
03492
03493 for (octave_idx_type k = 0; k < nc; k++)
03494 {
03495 if (work[k] != 0.)
03496 {
03497 octave_idx_type minr = nr;
03498 octave_idx_type mini = 0;
03499
03500 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
03501 if (perm[ridx(i)] < minr)
03502 {
03503 minr = perm[ridx(i)];
03504 mini = i;
03505 }
03506
03507 if (minr != k || data (mini) == 0.)
03508 {
03509 err = -2;
03510 goto triangular_error;
03511 }
03512
03513 Complex tmp = work[k] / data(mini);
03514 work[k] = tmp;
03515 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
03516 {
03517 if (i == mini)
03518 continue;
03519
03520 octave_idx_type iidx = perm[ridx(i)];
03521 work[iidx] = work[iidx] - tmp * data(i);
03522 }
03523 }
03524 }
03525
03526
03527
03528 octave_idx_type new_nnz = 0;
03529 for (octave_idx_type i = 0; i < nc; i++)
03530 if (work[i] != 0.)
03531 new_nnz++;
03532
03533 if (ii + new_nnz > x_nz)
03534 {
03535
03536 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
03537 retval.change_capacity (sz);
03538 x_nz = sz;
03539 }
03540
03541 for (octave_idx_type i = 0; i < nc; i++)
03542 if (work[i] != 0.)
03543 {
03544 retval.xridx(ii) = i;
03545 retval.xdata(ii++) = work[i];
03546 }
03547 retval.xcidx(j+1) = ii;
03548 }
03549
03550 retval.maybe_compress ();
03551
03552 if (calc_cond)
03553 {
03554
03555 for (octave_idx_type i = 0; i < nm; i++)
03556 work[i] = 0.;
03557
03558 for (octave_idx_type j = 0; j < nr; j++)
03559 {
03560 work[j] = 1.;
03561
03562 for (octave_idx_type k = 0; k < nc; k++)
03563 {
03564 if (work[k] != 0.)
03565 {
03566 octave_idx_type minr = nr;
03567 octave_idx_type mini = 0;
03568
03569 for (octave_idx_type i = cidx(k);
03570 i < cidx(k+1); i++)
03571 if (perm[ridx(i)] < minr)
03572 {
03573 minr = perm[ridx(i)];
03574 mini = i;
03575 }
03576
03577 Complex tmp = work[k] / data(mini);
03578 work[k] = tmp;
03579 for (octave_idx_type i = cidx(k);
03580 i < cidx(k+1); i++)
03581 {
03582 if (i == mini)
03583 continue;
03584
03585 octave_idx_type iidx = perm[ridx(i)];
03586 work[iidx] = work[iidx] - tmp * data(i);
03587 }
03588 }
03589 }
03590
03591 double atmp = 0;
03592 for (octave_idx_type i = j; i < nc; i++)
03593 {
03594 atmp += std::abs(work[i]);
03595 work[i] = 0.;
03596 }
03597 if (atmp > ainvnorm)
03598 ainvnorm = atmp;
03599 }
03600 rcond = 1. / ainvnorm / anorm;
03601 }
03602 }
03603 else
03604 {
03605 OCTAVE_LOCAL_BUFFER (Complex, work, nm);
03606
03607 for (octave_idx_type j = 0; j < b_nc; j++)
03608 {
03609 for (octave_idx_type i = 0; i < nm; i++)
03610 work[i] = 0.;
03611 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
03612 work[b.ridx(i)] = b.data(i);
03613
03614 for (octave_idx_type k = 0; k < nc; k++)
03615 {
03616 if (work[k] != 0.)
03617 {
03618 if (ridx(cidx(k)) != k ||
03619 data(cidx(k)) == 0.)
03620 {
03621 err = -2;
03622 goto triangular_error;
03623 }
03624
03625 Complex tmp = work[k] / data(cidx(k));
03626 work[k] = tmp;
03627 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
03628 {
03629 octave_idx_type iidx = ridx(i);
03630 work[iidx] = work[iidx] - tmp * data(i);
03631 }
03632 }
03633 }
03634
03635
03636
03637 octave_idx_type new_nnz = 0;
03638 for (octave_idx_type i = 0; i < nc; i++)
03639 if (work[i] != 0.)
03640 new_nnz++;
03641
03642 if (ii + new_nnz > x_nz)
03643 {
03644
03645 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
03646 retval.change_capacity (sz);
03647 x_nz = sz;
03648 }
03649
03650 for (octave_idx_type i = 0; i < nc; i++)
03651 if (work[i] != 0.)
03652 {
03653 retval.xridx(ii) = i;
03654 retval.xdata(ii++) = work[i];
03655 }
03656 retval.xcidx(j+1) = ii;
03657 }
03658
03659 retval.maybe_compress ();
03660
03661 if (calc_cond)
03662 {
03663
03664 for (octave_idx_type i = 0; i < nm; i++)
03665 work[i] = 0.;
03666
03667 for (octave_idx_type j = 0; j < nr; j++)
03668 {
03669 work[j] = 1.;
03670
03671 for (octave_idx_type k = j; k < nc; k++)
03672 {
03673
03674 if (work[k] != 0.)
03675 {
03676 Complex tmp = work[k] / data(cidx(k));
03677 work[k] = tmp;
03678 for (octave_idx_type i = cidx(k)+1;
03679 i < cidx(k+1); i++)
03680 {
03681 octave_idx_type iidx = ridx(i);
03682 work[iidx] = work[iidx] - tmp * data(i);
03683 }
03684 }
03685 }
03686 double atmp = 0;
03687 for (octave_idx_type i = j; i < nc; i++)
03688 {
03689 atmp += std::abs(work[i]);
03690 work[i] = 0.;
03691 }
03692 if (atmp > ainvnorm)
03693 ainvnorm = atmp;
03694 }
03695 rcond = 1. / ainvnorm / anorm;
03696 }
03697 }
03698
03699 triangular_error:
03700 if (err != 0)
03701 {
03702 if (sing_handler)
03703 {
03704 sing_handler (rcond);
03705 mattype.mark_as_rectangular ();
03706 }
03707 else
03708 (*current_liboctave_error_handler)
03709 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
03710 rcond);
03711 }
03712
03713 volatile double rcond_plus_one = rcond + 1.0;
03714
03715 if (rcond_plus_one == 1.0 || xisnan (rcond))
03716 {
03717 err = -2;
03718
03719 if (sing_handler)
03720 {
03721 sing_handler (rcond);
03722 mattype.mark_as_rectangular ();
03723 }
03724 else
03725 (*current_liboctave_error_handler)
03726 ("matrix singular to machine precision, rcond = %g",
03727 rcond);
03728 }
03729 }
03730 else
03731 (*current_liboctave_error_handler) ("incorrect matrix type");
03732 }
03733
03734 return retval;
03735 }
03736
03737 ComplexMatrix
03738 SparseComplexMatrix::trisolve (MatrixType &mattype, const Matrix& b,
03739 octave_idx_type& err, double& rcond,
03740 solve_singularity_handler sing_handler,
03741 bool calc_cond) const
03742 {
03743 ComplexMatrix retval;
03744
03745 octave_idx_type nr = rows ();
03746 octave_idx_type nc = cols ();
03747 err = 0;
03748
03749 if (nr != nc || nr != b.rows ())
03750 (*current_liboctave_error_handler)
03751 ("matrix dimension mismatch solution of linear equations");
03752 else if (nr == 0 || b.cols () == 0)
03753 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
03754 else if (calc_cond)
03755 (*current_liboctave_error_handler)
03756 ("calculation of condition number not implemented");
03757 else
03758 {
03759
03760 volatile int typ = mattype.type ();
03761 mattype.info ();
03762
03763 if (typ == MatrixType::Tridiagonal_Hermitian)
03764 {
03765 OCTAVE_LOCAL_BUFFER (double, D, nr);
03766 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
03767
03768 if (mattype.is_dense ())
03769 {
03770 octave_idx_type ii = 0;
03771
03772 for (octave_idx_type j = 0; j < nc-1; j++)
03773 {
03774 D[j] = std::real(data(ii++));
03775 DL[j] = data(ii);
03776 ii += 2;
03777 }
03778 D[nc-1] = std::real(data(ii));
03779 }
03780 else
03781 {
03782 D[0] = 0.;
03783 for (octave_idx_type i = 0; i < nr - 1; i++)
03784 {
03785 D[i+1] = 0.;
03786 DL[i] = 0.;
03787 }
03788
03789 for (octave_idx_type j = 0; j < nc; j++)
03790 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
03791 {
03792 if (ridx(i) == j)
03793 D[j] = std::real(data(i));
03794 else if (ridx(i) == j + 1)
03795 DL[j] = data(i);
03796 }
03797 }
03798
03799 octave_idx_type b_nc = b.cols();
03800 retval = ComplexMatrix (b);
03801 Complex *result = retval.fortran_vec ();
03802
03803 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
03804 b.rows(), err));
03805
03806 if (err != 0)
03807 {
03808 err = 0;
03809 mattype.mark_as_unsymmetric ();
03810 typ = MatrixType::Tridiagonal;
03811 }
03812 else
03813 rcond = 1.;
03814 }
03815
03816 if (typ == MatrixType::Tridiagonal)
03817 {
03818 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
03819 OCTAVE_LOCAL_BUFFER (Complex, D, nr);
03820 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
03821
03822 if (mattype.is_dense ())
03823 {
03824 octave_idx_type ii = 0;
03825
03826 for (octave_idx_type j = 0; j < nc-1; j++)
03827 {
03828 D[j] = data(ii++);
03829 DL[j] = data(ii++);
03830 DU[j] = data(ii++);
03831 }
03832 D[nc-1] = data(ii);
03833 }
03834 else
03835 {
03836 D[0] = 0.;
03837 for (octave_idx_type i = 0; i < nr - 1; i++)
03838 {
03839 D[i+1] = 0.;
03840 DL[i] = 0.;
03841 DU[i] = 0.;
03842 }
03843
03844 for (octave_idx_type j = 0; j < nc; j++)
03845 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
03846 {
03847 if (ridx(i) == j)
03848 D[j] = data(i);
03849 else if (ridx(i) == j + 1)
03850 DL[j] = data(i);
03851 else if (ridx(i) == j - 1)
03852 DU[j-1] = data(i);
03853 }
03854 }
03855
03856 octave_idx_type b_nc = b.cols();
03857 retval = ComplexMatrix (b);
03858 Complex *result = retval.fortran_vec ();
03859
03860 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
03861 b.rows(), err));
03862
03863 if (err != 0)
03864 {
03865 rcond = 0.;
03866 err = -2;
03867
03868 if (sing_handler)
03869 {
03870 sing_handler (rcond);
03871 mattype.mark_as_rectangular ();
03872 }
03873 else
03874 (*current_liboctave_error_handler)
03875 ("matrix singular to machine precision");
03876
03877 }
03878 else
03879 rcond = 1.;
03880 }
03881 else if (typ != MatrixType::Tridiagonal_Hermitian)
03882 (*current_liboctave_error_handler) ("incorrect matrix type");
03883 }
03884
03885 return retval;
03886 }
03887
03888 SparseComplexMatrix
03889 SparseComplexMatrix::trisolve (MatrixType &mattype, const SparseMatrix& b,
03890 octave_idx_type& err, double& rcond,
03891 solve_singularity_handler sing_handler,
03892 bool calc_cond) const
03893 {
03894 SparseComplexMatrix retval;
03895
03896 octave_idx_type nr = rows ();
03897 octave_idx_type nc = cols ();
03898 err = 0;
03899
03900 if (nr != nc || nr != b.rows ())
03901 (*current_liboctave_error_handler)
03902 ("matrix dimension mismatch solution of linear equations");
03903 else if (nr == 0 || b.cols () == 0)
03904 retval = SparseComplexMatrix (nc, b.cols ());
03905 else if (calc_cond)
03906 (*current_liboctave_error_handler)
03907 ("calculation of condition number not implemented");
03908 else
03909 {
03910
03911 int typ = mattype.type ();
03912 mattype.info ();
03913
03914
03915 if (typ == MatrixType::Tridiagonal ||
03916 typ == MatrixType::Tridiagonal_Hermitian)
03917 {
03918 OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
03919 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
03920 OCTAVE_LOCAL_BUFFER (Complex, D, nr);
03921 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
03922 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
03923 octave_idx_type *pipvt = ipvt.fortran_vec ();
03924
03925 if (mattype.is_dense ())
03926 {
03927 octave_idx_type ii = 0;
03928
03929 for (octave_idx_type j = 0; j < nc-1; j++)
03930 {
03931 D[j] = data(ii++);
03932 DL[j] = data(ii++);
03933 DU[j] = data(ii++);
03934 }
03935 D[nc-1] = data(ii);
03936 }
03937 else
03938 {
03939 D[0] = 0.;
03940 for (octave_idx_type i = 0; i < nr - 1; i++)
03941 {
03942 D[i+1] = 0.;
03943 DL[i] = 0.;
03944 DU[i] = 0.;
03945 }
03946
03947 for (octave_idx_type j = 0; j < nc; j++)
03948 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
03949 {
03950 if (ridx(i) == j)
03951 D[j] = data(i);
03952 else if (ridx(i) == j + 1)
03953 DL[j] = data(i);
03954 else if (ridx(i) == j - 1)
03955 DU[j-1] = data(i);
03956 }
03957 }
03958
03959 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
03960
03961 if (err != 0)
03962 {
03963 err = -2;
03964 rcond = 0.0;
03965
03966 if (sing_handler)
03967 {
03968 sing_handler (rcond);
03969 mattype.mark_as_rectangular ();
03970 }
03971 else
03972 (*current_liboctave_error_handler)
03973 ("matrix singular to machine precision");
03974
03975 }
03976 else
03977 {
03978 char job = 'N';
03979 volatile octave_idx_type x_nz = b.nnz ();
03980 octave_idx_type b_nc = b.cols ();
03981 retval = SparseComplexMatrix (nr, b_nc, x_nz);
03982 retval.xcidx(0) = 0;
03983 volatile octave_idx_type ii = 0;
03984 rcond = 1.0;
03985
03986 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
03987
03988 for (volatile octave_idx_type j = 0; j < b_nc; j++)
03989 {
03990 for (octave_idx_type i = 0; i < nr; i++)
03991 work[i] = 0.;
03992 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
03993 work[b.ridx(i)] = b.data(i);
03994
03995 F77_XFCN (zgttrs, ZGTTRS,
03996 (F77_CONST_CHAR_ARG2 (&job, 1),
03997 nr, 1, DL, D, DU, DU2, pipvt,
03998 work, b.rows (), err
03999 F77_CHAR_ARG_LEN (1)));
04000
04001
04002
04003 octave_idx_type new_nnz = 0;
04004 for (octave_idx_type i = 0; i < nr; i++)
04005 if (work[i] != 0.)
04006 new_nnz++;
04007
04008 if (ii + new_nnz > x_nz)
04009 {
04010
04011 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
04012 retval.change_capacity (sz);
04013 x_nz = sz;
04014 }
04015
04016 for (octave_idx_type i = 0; i < nr; i++)
04017 if (work[i] != 0.)
04018 {
04019 retval.xridx(ii) = i;
04020 retval.xdata(ii++) = work[i];
04021 }
04022 retval.xcidx(j+1) = ii;
04023 }
04024
04025 retval.maybe_compress ();
04026 }
04027 }
04028 else if (typ != MatrixType::Tridiagonal_Hermitian)
04029 (*current_liboctave_error_handler) ("incorrect matrix type");
04030 }
04031
04032 return retval;
04033 }
04034
04035 ComplexMatrix
04036 SparseComplexMatrix::trisolve (MatrixType &mattype, const ComplexMatrix& b,
04037 octave_idx_type& err, double& rcond,
04038 solve_singularity_handler sing_handler,
04039 bool calc_cond) const
04040 {
04041 ComplexMatrix retval;
04042
04043 octave_idx_type nr = rows ();
04044 octave_idx_type nc = cols ();
04045 err = 0;
04046
04047 if (nr != nc || nr != b.rows ())
04048 (*current_liboctave_error_handler)
04049 ("matrix dimension mismatch solution of linear equations");
04050 else if (nr == 0 || b.cols () == 0)
04051 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
04052 else if (calc_cond)
04053 (*current_liboctave_error_handler)
04054 ("calculation of condition number not implemented");
04055 else
04056 {
04057
04058 volatile int typ = mattype.type ();
04059 mattype.info ();
04060
04061 if (typ == MatrixType::Tridiagonal_Hermitian)
04062 {
04063 OCTAVE_LOCAL_BUFFER (double, D, nr);
04064 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
04065
04066 if (mattype.is_dense ())
04067 {
04068 octave_idx_type ii = 0;
04069
04070 for (octave_idx_type j = 0; j < nc-1; j++)
04071 {
04072 D[j] = std::real(data(ii++));
04073 DL[j] = data(ii);
04074 ii += 2;
04075 }
04076 D[nc-1] = std::real(data(ii));
04077 }
04078 else
04079 {
04080 D[0] = 0.;
04081 for (octave_idx_type i = 0; i < nr - 1; i++)
04082 {
04083 D[i+1] = 0.;
04084 DL[i] = 0.;
04085 }
04086
04087 for (octave_idx_type j = 0; j < nc; j++)
04088 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04089 {
04090 if (ridx(i) == j)
04091 D[j] = std::real (data(i));
04092 else if (ridx(i) == j + 1)
04093 DL[j] = data(i);
04094 }
04095 }
04096
04097 octave_idx_type b_nr = b.rows ();
04098 octave_idx_type b_nc = b.cols();
04099 rcond = 1.;
04100
04101 retval = ComplexMatrix (b);
04102 Complex *result = retval.fortran_vec ();
04103
04104 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
04105 b_nr, err));
04106
04107 if (err != 0)
04108 {
04109 err = 0;
04110 mattype.mark_as_unsymmetric ();
04111 typ = MatrixType::Tridiagonal;
04112 }
04113 }
04114
04115 if (typ == MatrixType::Tridiagonal)
04116 {
04117 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
04118 OCTAVE_LOCAL_BUFFER (Complex, D, nr);
04119 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
04120
04121 if (mattype.is_dense ())
04122 {
04123 octave_idx_type ii = 0;
04124
04125 for (octave_idx_type j = 0; j < nc-1; j++)
04126 {
04127 D[j] = data(ii++);
04128 DL[j] = data(ii++);
04129 DU[j] = data(ii++);
04130 }
04131 D[nc-1] = data(ii);
04132 }
04133 else
04134 {
04135 D[0] = 0.;
04136 for (octave_idx_type i = 0; i < nr - 1; i++)
04137 {
04138 D[i+1] = 0.;
04139 DL[i] = 0.;
04140 DU[i] = 0.;
04141 }
04142
04143 for (octave_idx_type j = 0; j < nc; j++)
04144 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04145 {
04146 if (ridx(i) == j)
04147 D[j] = data(i);
04148 else if (ridx(i) == j + 1)
04149 DL[j] = data(i);
04150 else if (ridx(i) == j - 1)
04151 DU[j-1] = data(i);
04152 }
04153 }
04154
04155 octave_idx_type b_nr = b.rows();
04156 octave_idx_type b_nc = b.cols();
04157 rcond = 1.;
04158
04159 retval = ComplexMatrix (b);
04160 Complex *result = retval.fortran_vec ();
04161
04162 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
04163 b_nr, err));
04164
04165 if (err != 0)
04166 {
04167 rcond = 0.;
04168 err = -2;
04169
04170 if (sing_handler)
04171 {
04172 sing_handler (rcond);
04173 mattype.mark_as_rectangular ();
04174 }
04175 else
04176 (*current_liboctave_error_handler)
04177 ("matrix singular to machine precision");
04178 }
04179 }
04180 else if (typ != MatrixType::Tridiagonal_Hermitian)
04181 (*current_liboctave_error_handler) ("incorrect matrix type");
04182 }
04183
04184 return retval;
04185 }
04186
04187 SparseComplexMatrix
04188 SparseComplexMatrix::trisolve (MatrixType &mattype,
04189 const SparseComplexMatrix& b,
04190 octave_idx_type& err, double& rcond,
04191 solve_singularity_handler sing_handler,
04192 bool calc_cond) const
04193 {
04194 SparseComplexMatrix retval;
04195
04196 octave_idx_type nr = rows ();
04197 octave_idx_type nc = cols ();
04198 err = 0;
04199
04200 if (nr != nc || nr != b.rows ())
04201 (*current_liboctave_error_handler)
04202 ("matrix dimension mismatch solution of linear equations");
04203 else if (nr == 0 || b.cols () == 0)
04204 retval = SparseComplexMatrix (nc, b.cols ());
04205 else if (calc_cond)
04206 (*current_liboctave_error_handler)
04207 ("calculation of condition number not implemented");
04208 else
04209 {
04210
04211 int typ = mattype.type ();
04212 mattype.info ();
04213
04214
04215 if (typ == MatrixType::Tridiagonal ||
04216 typ == MatrixType::Tridiagonal_Hermitian)
04217 {
04218 OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
04219 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
04220 OCTAVE_LOCAL_BUFFER (Complex, D, nr);
04221 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
04222 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
04223 octave_idx_type *pipvt = ipvt.fortran_vec ();
04224
04225 if (mattype.is_dense ())
04226 {
04227 octave_idx_type ii = 0;
04228
04229 for (octave_idx_type j = 0; j < nc-1; j++)
04230 {
04231 D[j] = data(ii++);
04232 DL[j] = data(ii++);
04233 DU[j] = data(ii++);
04234 }
04235 D[nc-1] = data(ii);
04236 }
04237 else
04238 {
04239 D[0] = 0.;
04240 for (octave_idx_type i = 0; i < nr - 1; i++)
04241 {
04242 D[i+1] = 0.;
04243 DL[i] = 0.;
04244 DU[i] = 0.;
04245 }
04246
04247 for (octave_idx_type j = 0; j < nc; j++)
04248 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04249 {
04250 if (ridx(i) == j)
04251 D[j] = data(i);
04252 else if (ridx(i) == j + 1)
04253 DL[j] = data(i);
04254 else if (ridx(i) == j - 1)
04255 DU[j-1] = data(i);
04256 }
04257 }
04258
04259 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
04260
04261 if (err != 0)
04262 {
04263 rcond = 0.0;
04264 err = -2;
04265
04266 if (sing_handler)
04267 {
04268 sing_handler (rcond);
04269 mattype.mark_as_rectangular ();
04270 }
04271 else
04272 (*current_liboctave_error_handler)
04273 ("matrix singular to machine precision");
04274 }
04275 else
04276 {
04277 rcond = 1.;
04278 char job = 'N';
04279 octave_idx_type b_nr = b.rows ();
04280 octave_idx_type b_nc = b.cols ();
04281 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
04282
04283
04284
04285 volatile octave_idx_type x_nz = b.nnz ();
04286 volatile octave_idx_type ii = 0;
04287 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
04288
04289 retval.xcidx(0) = 0;
04290 for (volatile octave_idx_type j = 0; j < b_nc; j++)
04291 {
04292
04293 for (octave_idx_type i = 0; i < b_nr; i++)
04294 Bx[i] = b (i,j);
04295
04296 F77_XFCN (zgttrs, ZGTTRS,
04297 (F77_CONST_CHAR_ARG2 (&job, 1),
04298 nr, 1, DL, D, DU, DU2, pipvt,
04299 Bx, b_nr, err
04300 F77_CHAR_ARG_LEN (1)));
04301
04302 if (err != 0)
04303 {
04304 (*current_liboctave_error_handler)
04305 ("SparseComplexMatrix::solve solve failed");
04306
04307 err = -1;
04308 break;
04309 }
04310
04311
04312
04313 octave_idx_type new_nnz = 0;
04314 for (octave_idx_type i = 0; i < nr; i++)
04315 if (Bx[i] != 0.)
04316 new_nnz++;
04317
04318 if (ii + new_nnz > x_nz)
04319 {
04320
04321 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
04322 retval.change_capacity (sz);
04323 x_nz = sz;
04324 }
04325
04326 for (octave_idx_type i = 0; i < nr; i++)
04327 if (Bx[i] != 0.)
04328 {
04329 retval.xridx(ii) = i;
04330 retval.xdata(ii++) = Bx[i];
04331 }
04332
04333 retval.xcidx(j+1) = ii;
04334 }
04335
04336 retval.maybe_compress ();
04337 }
04338 }
04339 else if (typ != MatrixType::Tridiagonal_Hermitian)
04340 (*current_liboctave_error_handler) ("incorrect matrix type");
04341 }
04342
04343 return retval;
04344 }
04345
04346 ComplexMatrix
04347 SparseComplexMatrix::bsolve (MatrixType &mattype, const Matrix& b,
04348 octave_idx_type& err, double& rcond,
04349 solve_singularity_handler sing_handler,
04350 bool calc_cond) const
04351 {
04352 ComplexMatrix retval;
04353
04354 octave_idx_type nr = rows ();
04355 octave_idx_type nc = cols ();
04356 err = 0;
04357
04358 if (nr != nc || nr != b.rows ())
04359 (*current_liboctave_error_handler)
04360 ("matrix dimension mismatch solution of linear equations");
04361 else if (nr == 0 || b.cols () == 0)
04362 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
04363 else
04364 {
04365
04366 volatile int typ = mattype.type ();
04367 mattype.info ();
04368
04369 if (typ == MatrixType::Banded_Hermitian)
04370 {
04371 octave_idx_type n_lower = mattype.nlower ();
04372 octave_idx_type ldm = n_lower + 1;
04373 ComplexMatrix m_band (ldm, nc);
04374 Complex *tmp_data = m_band.fortran_vec ();
04375
04376 if (! mattype.is_dense ())
04377 {
04378 octave_idx_type ii = 0;
04379
04380 for (octave_idx_type j = 0; j < ldm; j++)
04381 for (octave_idx_type i = 0; i < nc; i++)
04382 tmp_data[ii++] = 0.;
04383 }
04384
04385 for (octave_idx_type j = 0; j < nc; j++)
04386 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04387 {
04388 octave_idx_type ri = ridx (i);
04389 if (ri >= j)
04390 m_band(ri - j, j) = data(i);
04391 }
04392
04393
04394 double anorm;
04395 if (calc_cond)
04396 anorm = m_band.abs().sum().row(0).max();
04397
04398 char job = 'L';
04399 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
04400 nr, n_lower, tmp_data, ldm, err
04401 F77_CHAR_ARG_LEN (1)));
04402
04403 if (err != 0)
04404 {
04405 rcond = 0.0;
04406
04407
04408 mattype.mark_as_unsymmetric ();
04409 typ = MatrixType::Banded;
04410 err = 0;
04411 }
04412 else
04413 {
04414 if (calc_cond)
04415 {
04416 Array<Complex> z (dim_vector (2 * nr, 1));
04417 Complex *pz = z.fortran_vec ();
04418 Array<double> iz (dim_vector (nr, 1));
04419 double *piz = iz.fortran_vec ();
04420
04421 F77_XFCN (zpbcon, ZPBCON,
04422 (F77_CONST_CHAR_ARG2 (&job, 1),
04423 nr, n_lower, tmp_data, ldm,
04424 anorm, rcond, pz, piz, err
04425 F77_CHAR_ARG_LEN (1)));
04426
04427 if (err != 0)
04428 err = -2;
04429
04430 volatile double rcond_plus_one = rcond + 1.0;
04431
04432 if (rcond_plus_one == 1.0 || xisnan (rcond))
04433 {
04434 err = -2;
04435
04436 if (sing_handler)
04437 {
04438 sing_handler (rcond);
04439 mattype.mark_as_rectangular ();
04440 }
04441 else
04442 (*current_liboctave_error_handler)
04443 ("matrix singular to machine precision, rcond = %g",
04444 rcond);
04445 }
04446 }
04447 else
04448 rcond = 1.0;
04449
04450 if (err == 0)
04451 {
04452 retval = ComplexMatrix (b);
04453 Complex *result = retval.fortran_vec ();
04454
04455 octave_idx_type b_nc = b.cols ();
04456
04457 F77_XFCN (zpbtrs, ZPBTRS,
04458 (F77_CONST_CHAR_ARG2 (&job, 1),
04459 nr, n_lower, b_nc, tmp_data,
04460 ldm, result, b.rows(), err
04461 F77_CHAR_ARG_LEN (1)));
04462
04463 if (err != 0)
04464 {
04465 (*current_liboctave_error_handler)
04466 ("SparseMatrix::solve solve failed");
04467 err = -1;
04468 }
04469 }
04470 }
04471 }
04472
04473 if (typ == MatrixType::Banded)
04474 {
04475
04476 octave_idx_type n_upper = mattype.nupper ();
04477 octave_idx_type n_lower = mattype.nlower ();
04478 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
04479
04480 ComplexMatrix m_band (ldm, nc);
04481 Complex *tmp_data = m_band.fortran_vec ();
04482
04483 if (! mattype.is_dense ())
04484 {
04485 octave_idx_type ii = 0;
04486
04487 for (octave_idx_type j = 0; j < ldm; j++)
04488 for (octave_idx_type i = 0; i < nc; i++)
04489 tmp_data[ii++] = 0.;
04490 }
04491
04492 for (octave_idx_type j = 0; j < nc; j++)
04493 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04494 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
04495
04496
04497 double anorm;
04498 if (calc_cond)
04499 {
04500 for (octave_idx_type j = 0; j < nr; j++)
04501 {
04502 double atmp = 0.;
04503 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04504 atmp += std::abs(data(i));
04505 if (atmp > anorm)
04506 anorm = atmp;
04507 }
04508 }
04509
04510 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
04511 octave_idx_type *pipvt = ipvt.fortran_vec ();
04512
04513 F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data,
04514 ldm, pipvt, err));
04515
04516
04517
04518 if (err != 0)
04519 {
04520 rcond = 0.0;
04521 err = -2;
04522
04523 if (sing_handler)
04524 {
04525 sing_handler (rcond);
04526 mattype.mark_as_rectangular ();
04527 }
04528 else
04529 (*current_liboctave_error_handler)
04530 ("matrix singular to machine precision");
04531 }
04532 else
04533 {
04534 if (calc_cond)
04535 {
04536 char job = '1';
04537 Array<Complex> z (dim_vector (2 * nr, 1));
04538 Complex *pz = z.fortran_vec ();
04539 Array<double> iz (dim_vector (nr, 1));
04540 double *piz = iz.fortran_vec ();
04541
04542 F77_XFCN (zgbcon, ZGBCON,
04543 (F77_CONST_CHAR_ARG2 (&job, 1),
04544 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
04545 anorm, rcond, pz, piz, err
04546 F77_CHAR_ARG_LEN (1)));
04547
04548 if (err != 0)
04549 err = -2;
04550
04551 volatile double rcond_plus_one = rcond + 1.0;
04552
04553 if (rcond_plus_one == 1.0 || xisnan (rcond))
04554 {
04555 err = -2;
04556
04557 if (sing_handler)
04558 {
04559 sing_handler (rcond);
04560 mattype.mark_as_rectangular ();
04561 }
04562 else
04563 (*current_liboctave_error_handler)
04564 ("matrix singular to machine precision, rcond = %g",
04565 rcond);
04566 }
04567 }
04568 else
04569 rcond = 1.;
04570
04571 if (err == 0)
04572 {
04573 retval = ComplexMatrix (b);
04574 Complex *result = retval.fortran_vec ();
04575
04576 octave_idx_type b_nc = b.cols ();
04577
04578 char job = 'N';
04579 F77_XFCN (zgbtrs, ZGBTRS,
04580 (F77_CONST_CHAR_ARG2 (&job, 1),
04581 nr, n_lower, n_upper, b_nc, tmp_data,
04582 ldm, pipvt, result, b.rows(), err
04583 F77_CHAR_ARG_LEN (1)));
04584 }
04585 }
04586 }
04587 else if (typ != MatrixType::Banded_Hermitian)
04588 (*current_liboctave_error_handler) ("incorrect matrix type");
04589 }
04590
04591 return retval;
04592 }
04593
04594 SparseComplexMatrix
04595 SparseComplexMatrix::bsolve (MatrixType &mattype, const SparseMatrix& b,
04596 octave_idx_type& err, double& rcond,
04597 solve_singularity_handler sing_handler,
04598 bool calc_cond) const
04599 {
04600 SparseComplexMatrix retval;
04601
04602 octave_idx_type nr = rows ();
04603 octave_idx_type nc = cols ();
04604 err = 0;
04605
04606 if (nr != nc || nr != b.rows ())
04607 (*current_liboctave_error_handler)
04608 ("matrix dimension mismatch solution of linear equations");
04609 else if (nr == 0 || b.cols () == 0)
04610 retval = SparseComplexMatrix (nc, b.cols ());
04611 else
04612 {
04613
04614 volatile int typ = mattype.type ();
04615 mattype.info ();
04616
04617 if (typ == MatrixType::Banded_Hermitian)
04618 {
04619 octave_idx_type n_lower = mattype.nlower ();
04620 octave_idx_type ldm = n_lower + 1;
04621
04622 ComplexMatrix m_band (ldm, nc);
04623 Complex *tmp_data = m_band.fortran_vec ();
04624
04625 if (! mattype.is_dense ())
04626 {
04627 octave_idx_type ii = 0;
04628
04629 for (octave_idx_type j = 0; j < ldm; j++)
04630 for (octave_idx_type i = 0; i < nc; i++)
04631 tmp_data[ii++] = 0.;
04632 }
04633
04634 for (octave_idx_type j = 0; j < nc; j++)
04635 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04636 {
04637 octave_idx_type ri = ridx (i);
04638 if (ri >= j)
04639 m_band(ri - j, j) = data(i);
04640 }
04641
04642
04643 double anorm;
04644 if (calc_cond)
04645 anorm = m_band.abs().sum().row(0).max();
04646
04647 char job = 'L';
04648 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
04649 nr, n_lower, tmp_data, ldm, err
04650 F77_CHAR_ARG_LEN (1)));
04651
04652 if (err != 0)
04653 {
04654 rcond = 0.0;
04655 mattype.mark_as_unsymmetric ();
04656 typ = MatrixType::Banded;
04657 err = 0;
04658 }
04659 else
04660 {
04661 if (calc_cond)
04662 {
04663 Array<Complex> z (dim_vector (2 * nr, 1));
04664 Complex *pz = z.fortran_vec ();
04665 Array<double> iz (dim_vector (nr, 1));
04666 double *piz = iz.fortran_vec ();
04667
04668 F77_XFCN (zpbcon, ZPBCON,
04669 (F77_CONST_CHAR_ARG2 (&job, 1),
04670 nr, n_lower, tmp_data, ldm,
04671 anorm, rcond, pz, piz, err
04672 F77_CHAR_ARG_LEN (1)));
04673
04674 if (err != 0)
04675 err = -2;
04676
04677 volatile double rcond_plus_one = rcond + 1.0;
04678
04679 if (rcond_plus_one == 1.0 || xisnan (rcond))
04680 {
04681 err = -2;
04682
04683 if (sing_handler)
04684 {
04685 sing_handler (rcond);
04686 mattype.mark_as_rectangular ();
04687 }
04688 else
04689 (*current_liboctave_error_handler)
04690 ("matrix singular to machine precision, rcond = %g",
04691 rcond);
04692 }
04693 }
04694 else
04695 rcond = 1.0;
04696
04697 if (err == 0)
04698 {
04699 octave_idx_type b_nr = b.rows ();
04700 octave_idx_type b_nc = b.cols ();
04701 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
04702
04703
04704
04705 volatile octave_idx_type x_nz = b.nnz ();
04706 volatile octave_idx_type ii = 0;
04707 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
04708
04709 retval.xcidx(0) = 0;
04710 for (volatile octave_idx_type j = 0; j < b_nc; j++)
04711 {
04712 for (octave_idx_type i = 0; i < b_nr; i++)
04713 Bx[i] = b.elem (i, j);
04714
04715 F77_XFCN (zpbtrs, ZPBTRS,
04716 (F77_CONST_CHAR_ARG2 (&job, 1),
04717 nr, n_lower, 1, tmp_data,
04718 ldm, Bx, b_nr, err
04719 F77_CHAR_ARG_LEN (1)));
04720
04721 if (err != 0)
04722 {
04723 (*current_liboctave_error_handler)
04724 ("SparseComplexMatrix::solve solve failed");
04725 err = -1;
04726 break;
04727 }
04728
04729 for (octave_idx_type i = 0; i < b_nr; i++)
04730 {
04731 Complex tmp = Bx[i];
04732 if (tmp != 0.0)
04733 {
04734 if (ii == x_nz)
04735 {
04736
04737 octave_idx_type sz = x_nz *
04738 (b_nc - j) / b_nc;
04739 sz = (sz > 10 ? sz : 10) + x_nz;
04740 retval.change_capacity (sz);
04741 x_nz = sz;
04742 }
04743 retval.xdata(ii) = tmp;
04744 retval.xridx(ii++) = i;
04745 }
04746 }
04747 retval.xcidx(j+1) = ii;
04748 }
04749
04750 retval.maybe_compress ();
04751 }
04752 }
04753 }
04754
04755 if (typ == MatrixType::Banded)
04756 {
04757
04758 octave_idx_type n_upper = mattype.nupper ();
04759 octave_idx_type n_lower = mattype.nlower ();
04760 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
04761
04762 ComplexMatrix m_band (ldm, nc);
04763 Complex *tmp_data = m_band.fortran_vec ();
04764
04765 if (! mattype.is_dense ())
04766 {
04767 octave_idx_type ii = 0;
04768
04769 for (octave_idx_type j = 0; j < ldm; j++)
04770 for (octave_idx_type i = 0; i < nc; i++)
04771 tmp_data[ii++] = 0.;
04772 }
04773
04774 for (octave_idx_type j = 0; j < nc; j++)
04775 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04776 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
04777
04778
04779 double anorm;
04780 if (calc_cond)
04781 {
04782 for (octave_idx_type j = 0; j < nr; j++)
04783 {
04784 double atmp = 0.;
04785 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04786 atmp += std::abs(data(i));
04787 if (atmp > anorm)
04788 anorm = atmp;
04789 }
04790 }
04791
04792 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
04793 octave_idx_type *pipvt = ipvt.fortran_vec ();
04794
04795 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
04796 ldm, pipvt, err));
04797
04798 if (err != 0)
04799 {
04800 rcond = 0.0;
04801 err = -2;
04802
04803 if (sing_handler)
04804 {
04805 sing_handler (rcond);
04806 mattype.mark_as_rectangular ();
04807 }
04808 else
04809 (*current_liboctave_error_handler)
04810 ("matrix singular to machine precision");
04811
04812 }
04813 else
04814 {
04815 if (calc_cond)
04816 {
04817 char job = '1';
04818 Array<Complex> z (dim_vector (2 * nr, 1));
04819 Complex *pz = z.fortran_vec ();
04820 Array<double> iz (dim_vector (nr, 1));
04821 double *piz = iz.fortran_vec ();
04822
04823 F77_XFCN (zgbcon, ZGBCON,
04824 (F77_CONST_CHAR_ARG2 (&job, 1),
04825 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
04826 anorm, rcond, pz, piz, err
04827 F77_CHAR_ARG_LEN (1)));
04828
04829 if (err != 0)
04830 err = -2;
04831
04832 volatile double rcond_plus_one = rcond + 1.0;
04833
04834 if (rcond_plus_one == 1.0 || xisnan (rcond))
04835 {
04836 err = -2;
04837
04838 if (sing_handler)
04839 {
04840 sing_handler (rcond);
04841 mattype.mark_as_rectangular ();
04842 }
04843 else
04844 (*current_liboctave_error_handler)
04845 ("matrix singular to machine precision, rcond = %g",
04846 rcond);
04847 }
04848 }
04849 else
04850 rcond = 1.;
04851
04852 if (err == 0)
04853 {
04854 char job = 'N';
04855 volatile octave_idx_type x_nz = b.nnz ();
04856 octave_idx_type b_nc = b.cols ();
04857 retval = SparseComplexMatrix (nr, b_nc, x_nz);
04858 retval.xcidx(0) = 0;
04859 volatile octave_idx_type ii = 0;
04860
04861 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
04862
04863 for (volatile octave_idx_type j = 0; j < b_nc; j++)
04864 {
04865 for (octave_idx_type i = 0; i < nr; i++)
04866 work[i] = 0.;
04867 for (octave_idx_type i = b.cidx(j);
04868 i < b.cidx(j+1); i++)
04869 work[b.ridx(i)] = b.data(i);
04870
04871 F77_XFCN (zgbtrs, ZGBTRS,
04872 (F77_CONST_CHAR_ARG2 (&job, 1),
04873 nr, n_lower, n_upper, 1, tmp_data,
04874 ldm, pipvt, work, b.rows (), err
04875 F77_CHAR_ARG_LEN (1)));
04876
04877
04878
04879 octave_idx_type new_nnz = 0;
04880 for (octave_idx_type i = 0; i < nr; i++)
04881 if (work[i] != 0.)
04882 new_nnz++;
04883
04884 if (ii + new_nnz > x_nz)
04885 {
04886
04887 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
04888 retval.change_capacity (sz);
04889 x_nz = sz;
04890 }
04891
04892 for (octave_idx_type i = 0; i < nr; i++)
04893 if (work[i] != 0.)
04894 {
04895 retval.xridx(ii) = i;
04896 retval.xdata(ii++) = work[i];
04897 }
04898 retval.xcidx(j+1) = ii;
04899 }
04900
04901 retval.maybe_compress ();
04902 }
04903 }
04904 }
04905 else if (typ != MatrixType::Banded_Hermitian)
04906 (*current_liboctave_error_handler) ("incorrect matrix type");
04907 }
04908
04909 return retval;
04910 }
04911
04912 ComplexMatrix
04913 SparseComplexMatrix::bsolve (MatrixType &mattype, const ComplexMatrix& b,
04914 octave_idx_type& err, double& rcond,
04915 solve_singularity_handler sing_handler,
04916 bool calc_cond) const
04917 {
04918 ComplexMatrix retval;
04919
04920 octave_idx_type nr = rows ();
04921 octave_idx_type nc = cols ();
04922 err = 0;
04923
04924 if (nr != nc || nr != b.rows ())
04925 (*current_liboctave_error_handler)
04926 ("matrix dimension mismatch solution of linear equations");
04927 else if (nr == 0 || b.cols () == 0)
04928 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
04929 else
04930 {
04931
04932 volatile int typ = mattype.type ();
04933 mattype.info ();
04934
04935 if (typ == MatrixType::Banded_Hermitian)
04936 {
04937 octave_idx_type n_lower = mattype.nlower ();
04938 octave_idx_type ldm = n_lower + 1;
04939
04940 ComplexMatrix m_band (ldm, nc);
04941 Complex *tmp_data = m_band.fortran_vec ();
04942
04943 if (! mattype.is_dense ())
04944 {
04945 octave_idx_type ii = 0;
04946
04947 for (octave_idx_type j = 0; j < ldm; j++)
04948 for (octave_idx_type i = 0; i < nc; i++)
04949 tmp_data[ii++] = 0.;
04950 }
04951
04952 for (octave_idx_type j = 0; j < nc; j++)
04953 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
04954 {
04955 octave_idx_type ri = ridx (i);
04956 if (ri >= j)
04957 m_band(ri - j, j) = data(i);
04958 }
04959
04960
04961 double anorm;
04962 if (calc_cond)
04963 anorm = m_band.abs().sum().row(0).max();
04964
04965 char job = 'L';
04966 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
04967 nr, n_lower, tmp_data, ldm, err
04968 F77_CHAR_ARG_LEN (1)));
04969
04970 if (err != 0)
04971 {
04972
04973
04974 rcond = 0.0;
04975 mattype.mark_as_unsymmetric ();
04976 typ = MatrixType::Banded;
04977 err = 0;
04978 }
04979 else
04980 {
04981 if (calc_cond)
04982 {
04983 Array<Complex> z (dim_vector (2 * nr, 1));
04984 Complex *pz = z.fortran_vec ();
04985 Array<double> iz (dim_vector (nr, 1));
04986 double *piz = iz.fortran_vec ();
04987
04988 F77_XFCN (zpbcon, ZPBCON,
04989 (F77_CONST_CHAR_ARG2 (&job, 1),
04990 nr, n_lower, tmp_data, ldm,
04991 anorm, rcond, pz, piz, err
04992 F77_CHAR_ARG_LEN (1)));
04993
04994 if (err != 0)
04995 err = -2;
04996
04997 volatile double rcond_plus_one = rcond + 1.0;
04998
04999 if (rcond_plus_one == 1.0 || xisnan (rcond))
05000 {
05001 err = -2;
05002
05003 if (sing_handler)
05004 {
05005 sing_handler (rcond);
05006 mattype.mark_as_rectangular ();
05007 }
05008 else
05009 (*current_liboctave_error_handler)
05010 ("matrix singular to machine precision, rcond = %g",
05011 rcond);
05012 }
05013 }
05014 else
05015 rcond = 1.0;
05016
05017 if (err == 0)
05018 {
05019 octave_idx_type b_nr = b.rows ();
05020 octave_idx_type b_nc = b.cols ();
05021 retval = ComplexMatrix (b);
05022 Complex *result = retval.fortran_vec ();
05023
05024 F77_XFCN (zpbtrs, ZPBTRS,
05025 (F77_CONST_CHAR_ARG2 (&job, 1),
05026 nr, n_lower, b_nc, tmp_data,
05027 ldm, result, b_nr, err
05028 F77_CHAR_ARG_LEN (1)));
05029
05030 if (err != 0)
05031 {
05032 (*current_liboctave_error_handler)
05033 ("SparseComplexMatrix::solve solve failed");
05034 err = -1;
05035 }
05036 }
05037 }
05038 }
05039
05040 if (typ == MatrixType::Banded)
05041 {
05042
05043 octave_idx_type n_upper = mattype.nupper ();
05044 octave_idx_type n_lower = mattype.nlower ();
05045 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
05046
05047 ComplexMatrix m_band (ldm, nc);
05048 Complex *tmp_data = m_band.fortran_vec ();
05049
05050 if (! mattype.is_dense ())
05051 {
05052 octave_idx_type ii = 0;
05053
05054 for (octave_idx_type j = 0; j < ldm; j++)
05055 for (octave_idx_type i = 0; i < nc; i++)
05056 tmp_data[ii++] = 0.;
05057 }
05058
05059 for (octave_idx_type j = 0; j < nc; j++)
05060 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
05061 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
05062
05063
05064 double anorm;
05065 if (calc_cond)
05066 {
05067 for (octave_idx_type j = 0; j < nr; j++)
05068 {
05069 double atmp = 0.;
05070 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
05071 atmp += std::abs(data(i));
05072 if (atmp > anorm)
05073 anorm = atmp;
05074 }
05075 }
05076
05077 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
05078 octave_idx_type *pipvt = ipvt.fortran_vec ();
05079
05080 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
05081 ldm, pipvt, err));
05082
05083 if (err != 0)
05084 {
05085 err = -2;
05086 rcond = 0.0;
05087
05088 if (sing_handler)
05089 {
05090 sing_handler (rcond);
05091 mattype.mark_as_rectangular ();
05092 }
05093 else
05094 (*current_liboctave_error_handler)
05095 ("matrix singular to machine precision");
05096 }
05097 else
05098 {
05099 if (calc_cond)
05100 {
05101 char job = '1';
05102 Array<Complex> z (dim_vector (2 * nr, 1));
05103 Complex *pz = z.fortran_vec ();
05104 Array<double> iz (dim_vector (nr, 1));
05105 double *piz = iz.fortran_vec ();
05106
05107 F77_XFCN (zgbcon, ZGBCON,
05108 (F77_CONST_CHAR_ARG2 (&job, 1),
05109 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
05110 anorm, rcond, pz, piz, err
05111 F77_CHAR_ARG_LEN (1)));
05112
05113 if (err != 0)
05114 err = -2;
05115
05116 volatile double rcond_plus_one = rcond + 1.0;
05117
05118 if (rcond_plus_one == 1.0 || xisnan (rcond))
05119 {
05120 err = -2;
05121
05122 if (sing_handler)
05123 {
05124 sing_handler (rcond);
05125 mattype.mark_as_rectangular ();
05126 }
05127 else
05128 (*current_liboctave_error_handler)
05129 ("matrix singular to machine precision, rcond = %g",
05130 rcond);
05131 }
05132 }
05133 else
05134 rcond = 1.;
05135
05136 if (err == 0)
05137 {
05138 char job = 'N';
05139 octave_idx_type b_nc = b.cols ();
05140 retval = ComplexMatrix (b);
05141 Complex *result = retval.fortran_vec ();
05142
05143 F77_XFCN (zgbtrs, ZGBTRS,
05144 (F77_CONST_CHAR_ARG2 (&job, 1),
05145 nr, n_lower, n_upper, b_nc, tmp_data,
05146 ldm, pipvt, result, b.rows (), err
05147 F77_CHAR_ARG_LEN (1)));
05148 }
05149 }
05150 }
05151 else if (typ != MatrixType::Banded_Hermitian)
05152 (*current_liboctave_error_handler) ("incorrect matrix type");
05153 }
05154
05155 return retval;
05156 }
05157
05158 SparseComplexMatrix
05159 SparseComplexMatrix::bsolve (MatrixType &mattype, const SparseComplexMatrix& b,
05160 octave_idx_type& err, double& rcond,
05161 solve_singularity_handler sing_handler,
05162 bool calc_cond) const
05163 {
05164 SparseComplexMatrix retval;
05165
05166 octave_idx_type nr = rows ();
05167 octave_idx_type nc = cols ();
05168 err = 0;
05169
05170 if (nr != nc || nr != b.rows ())
05171 (*current_liboctave_error_handler)
05172 ("matrix dimension mismatch solution of linear equations");
05173 else if (nr == 0 || b.cols () == 0)
05174 retval = SparseComplexMatrix (nc, b.cols ());
05175 else
05176 {
05177
05178 volatile int typ = mattype.type ();
05179 mattype.info ();
05180
05181 if (typ == MatrixType::Banded_Hermitian)
05182 {
05183 octave_idx_type n_lower = mattype.nlower ();
05184 octave_idx_type ldm = n_lower + 1;
05185
05186 ComplexMatrix m_band (ldm, nc);
05187 Complex *tmp_data = m_band.fortran_vec ();
05188
05189 if (! mattype.is_dense ())
05190 {
05191 octave_idx_type ii = 0;
05192
05193 for (octave_idx_type j = 0; j < ldm; j++)
05194 for (octave_idx_type i = 0; i < nc; i++)
05195 tmp_data[ii++] = 0.;
05196 }
05197
05198 for (octave_idx_type j = 0; j < nc; j++)
05199 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
05200 {
05201 octave_idx_type ri = ridx (i);
05202 if (ri >= j)
05203 m_band(ri - j, j) = data(i);
05204 }
05205
05206
05207 double anorm;
05208 if (calc_cond)
05209 anorm = m_band.abs().sum().row(0).max();
05210
05211 char job = 'L';
05212 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
05213 nr, n_lower, tmp_data, ldm, err
05214 F77_CHAR_ARG_LEN (1)));
05215
05216 if (err != 0)
05217 {
05218
05219
05220 mattype.mark_as_unsymmetric ();
05221 typ = MatrixType::Banded;
05222
05223 rcond = 0.0;
05224 err = 0;
05225 }
05226 else
05227 {
05228 if (calc_cond)
05229 {
05230 Array<Complex> z (dim_vector (2 * nr, 1));
05231 Complex *pz = z.fortran_vec ();
05232 Array<double> iz (dim_vector (nr, 1));
05233 double *piz = iz.fortran_vec ();
05234
05235 F77_XFCN (zpbcon, ZPBCON,
05236 (F77_CONST_CHAR_ARG2 (&job, 1),
05237 nr, n_lower, tmp_data, ldm,
05238 anorm, rcond, pz, piz, err
05239 F77_CHAR_ARG_LEN (1)));
05240
05241 if (err != 0)
05242 err = -2;
05243
05244 volatile double rcond_plus_one = rcond + 1.0;
05245
05246 if (rcond_plus_one == 1.0 || xisnan (rcond))
05247 {
05248 err = -2;
05249
05250 if (sing_handler)
05251 {
05252 sing_handler (rcond);
05253 mattype.mark_as_rectangular ();
05254 }
05255 else
05256 (*current_liboctave_error_handler)
05257 ("matrix singular to machine precision, rcond = %g",
05258 rcond);
05259 }
05260 }
05261 else
05262 rcond = 1.0;
05263
05264 if (err == 0)
05265 {
05266 octave_idx_type b_nr = b.rows ();
05267 octave_idx_type b_nc = b.cols ();
05268 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
05269
05270
05271
05272 volatile octave_idx_type x_nz = b.nnz ();
05273 volatile octave_idx_type ii = 0;
05274 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
05275
05276 retval.xcidx(0) = 0;
05277 for (volatile octave_idx_type j = 0; j < b_nc; j++)
05278 {
05279
05280 for (octave_idx_type i = 0; i < b_nr; i++)
05281 Bx[i] = b (i,j);
05282
05283 F77_XFCN (zpbtrs, ZPBTRS,
05284 (F77_CONST_CHAR_ARG2 (&job, 1),
05285 nr, n_lower, 1, tmp_data,
05286 ldm, Bx, b_nr, err
05287 F77_CHAR_ARG_LEN (1)));
05288
05289 if (err != 0)
05290 {
05291 (*current_liboctave_error_handler)
05292 ("SparseMatrix::solve solve failed");
05293 err = -1;
05294 break;
05295 }
05296
05297
05298
05299 octave_idx_type new_nnz = 0;
05300 for (octave_idx_type i = 0; i < nr; i++)
05301 if (Bx[i] != 0.)
05302 new_nnz++;
05303
05304 if (ii + new_nnz > x_nz)
05305 {
05306
05307 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
05308 retval.change_capacity (sz);
05309 x_nz = sz;
05310 }
05311
05312 for (octave_idx_type i = 0; i < nr; i++)
05313 if (Bx[i] != 0.)
05314 {
05315 retval.xridx(ii) = i;
05316 retval.xdata(ii++) = Bx[i];
05317 }
05318
05319 retval.xcidx(j+1) = ii;
05320 }
05321
05322 retval.maybe_compress ();
05323 }
05324 }
05325 }
05326
05327 if (typ == MatrixType::Banded)
05328 {
05329
05330 octave_idx_type n_upper = mattype.nupper ();
05331 octave_idx_type n_lower = mattype.nlower ();
05332 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
05333
05334 ComplexMatrix m_band (ldm, nc);
05335 Complex *tmp_data = m_band.fortran_vec ();
05336
05337 if (! mattype.is_dense ())
05338 {
05339 octave_idx_type ii = 0;
05340
05341 for (octave_idx_type j = 0; j < ldm; j++)
05342 for (octave_idx_type i = 0; i < nc; i++)
05343 tmp_data[ii++] = 0.;
05344 }
05345
05346 for (octave_idx_type j = 0; j < nc; j++)
05347 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
05348 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
05349
05350
05351 double anorm;
05352 if (calc_cond)
05353 {
05354 for (octave_idx_type j = 0; j < nr; j++)
05355 {
05356 double atmp = 0.;
05357 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
05358 atmp += std::abs(data(i));
05359 if (atmp > anorm)
05360 anorm = atmp;
05361 }
05362 }
05363
05364 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
05365 octave_idx_type *pipvt = ipvt.fortran_vec ();
05366
05367 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
05368 ldm, pipvt, err));
05369
05370 if (err != 0)
05371 {
05372 err = -2;
05373 rcond = 0.0;
05374
05375 if (sing_handler)
05376 {
05377 sing_handler (rcond);
05378 mattype.mark_as_rectangular ();
05379 }
05380 else
05381 (*current_liboctave_error_handler)
05382 ("matrix singular to machine precision");
05383
05384 }
05385 else
05386 {
05387 if (calc_cond)
05388 {
05389 char job = '1';
05390 Array<Complex> z (dim_vector (2 * nr, 1));
05391 Complex *pz = z.fortran_vec ();
05392 Array<double> iz (dim_vector (nr, 1));
05393 double *piz = iz.fortran_vec ();
05394
05395 F77_XFCN (zgbcon, ZGBCON,
05396 (F77_CONST_CHAR_ARG2 (&job, 1),
05397 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
05398 anorm, rcond, pz, piz, err
05399 F77_CHAR_ARG_LEN (1)));
05400
05401 if (err != 0)
05402 err = -2;
05403
05404 volatile double rcond_plus_one = rcond + 1.0;
05405
05406 if (rcond_plus_one == 1.0 || xisnan (rcond))
05407 {
05408 err = -2;
05409
05410 if (sing_handler)
05411 {
05412 sing_handler (rcond);
05413 mattype.mark_as_rectangular ();
05414 }
05415 else
05416 (*current_liboctave_error_handler)
05417 ("matrix singular to machine precision, rcond = %g",
05418 rcond);
05419 }
05420 }
05421 else
05422 rcond = 1.;
05423
05424 if (err == 0)
05425 {
05426 char job = 'N';
05427 volatile octave_idx_type x_nz = b.nnz ();
05428 octave_idx_type b_nc = b.cols ();
05429 retval = SparseComplexMatrix (nr, b_nc, x_nz);
05430 retval.xcidx(0) = 0;
05431 volatile octave_idx_type ii = 0;
05432
05433 OCTAVE_LOCAL_BUFFER (Complex, Bx, nr);
05434
05435 for (volatile octave_idx_type j = 0; j < b_nc; j++)
05436 {
05437 for (octave_idx_type i = 0; i < nr; i++)
05438 Bx[i] = 0.;
05439
05440 for (octave_idx_type i = b.cidx(j);
05441 i < b.cidx(j+1); i++)
05442 Bx[b.ridx(i)] = b.data(i);
05443
05444 F77_XFCN (zgbtrs, ZGBTRS,
05445 (F77_CONST_CHAR_ARG2 (&job, 1),
05446 nr, n_lower, n_upper, 1, tmp_data,
05447 ldm, pipvt, Bx, b.rows (), err
05448 F77_CHAR_ARG_LEN (1)));
05449
05450
05451
05452 octave_idx_type new_nnz = 0;
05453 for (octave_idx_type i = 0; i < nr; i++)
05454 if (Bx[i] != 0.)
05455 new_nnz++;
05456
05457 if (ii + new_nnz > x_nz)
05458 {
05459
05460 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
05461 retval.change_capacity (sz);
05462 x_nz = sz;
05463 }
05464
05465 for (octave_idx_type i = 0; i < nr; i++)
05466 if (Bx[i] != 0.)
05467 {
05468 retval.xridx(ii) = i;
05469 retval.xdata(ii++) = Bx[i];
05470 }
05471 retval.xcidx(j+1) = ii;
05472 }
05473
05474 retval.maybe_compress ();
05475 }
05476 }
05477 }
05478 else if (typ != MatrixType::Banded_Hermitian)
05479 (*current_liboctave_error_handler) ("incorrect matrix type");
05480 }
05481
05482 return retval;
05483 }
05484
05485 void *
05486 SparseComplexMatrix::factorize (octave_idx_type& err, double &rcond,
05487 Matrix &Control, Matrix &Info,
05488 solve_singularity_handler sing_handler,
05489 bool calc_cond) const
05490 {
05491
05492 void *Numeric = 0;
05493 err = 0;
05494
05495 #ifdef HAVE_UMFPACK
05496
05497 Control = Matrix (UMFPACK_CONTROL, 1);
05498 double *control = Control.fortran_vec ();
05499 UMFPACK_ZNAME (defaults) (control);
05500
05501 double tmp = octave_sparse_params::get_key ("spumoni");
05502 if (!xisnan (tmp))
05503 Control (UMFPACK_PRL) = tmp;
05504 tmp = octave_sparse_params::get_key ("piv_tol");
05505 if (!xisnan (tmp))
05506 {
05507 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
05508 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
05509 }
05510
05511
05512 tmp = octave_sparse_params::get_key ("autoamd");
05513 if (!xisnan (tmp))
05514 Control (UMFPACK_FIXQ) = tmp;
05515
05516 UMFPACK_ZNAME (report_control) (control);
05517
05518 const octave_idx_type *Ap = cidx ();
05519 const octave_idx_type *Ai = ridx ();
05520 const Complex *Ax = data ();
05521 octave_idx_type nr = rows ();
05522 octave_idx_type nc = cols ();
05523
05524 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
05525 reinterpret_cast<const double *> (Ax),
05526 0, 1, control);
05527
05528 void *Symbolic;
05529 Info = Matrix (1, UMFPACK_INFO);
05530 double *info = Info.fortran_vec ();
05531 int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
05532 reinterpret_cast<const double *> (Ax),
05533 0, 0, &Symbolic, control, info);
05534
05535 if (status < 0)
05536 {
05537 (*current_liboctave_error_handler)
05538 ("SparseComplexMatrix::solve symbolic factorization failed");
05539 err = -1;
05540
05541 UMFPACK_ZNAME (report_status) (control, status);
05542 UMFPACK_ZNAME (report_info) (control, info);
05543
05544 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
05545 }
05546 else
05547 {
05548 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
05549
05550 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
05551 reinterpret_cast<const double *> (Ax), 0,
05552 Symbolic, &Numeric, control, info) ;
05553 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
05554
05555 if (calc_cond)
05556 rcond = Info (UMFPACK_RCOND);
05557 else
05558 rcond = 1.;
05559 volatile double rcond_plus_one = rcond + 1.0;
05560
05561 if (status == UMFPACK_WARNING_singular_matrix ||
05562 rcond_plus_one == 1.0 || xisnan (rcond))
05563 {
05564 UMFPACK_ZNAME (report_numeric) (Numeric, control);
05565
05566 err = -2;
05567
05568 if (sing_handler)
05569 sing_handler (rcond);
05570 else
05571 (*current_liboctave_error_handler)
05572 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
05573 rcond);
05574
05575 }
05576 else if (status < 0)
05577 {
05578 (*current_liboctave_error_handler)
05579 ("SparseComplexMatrix::solve numeric factorization failed");
05580
05581 UMFPACK_ZNAME (report_status) (control, status);
05582 UMFPACK_ZNAME (report_info) (control, info);
05583
05584 err = -1;
05585 }
05586 else
05587 {
05588 UMFPACK_ZNAME (report_numeric) (Numeric, control);
05589 }
05590 }
05591
05592 if (err != 0)
05593 UMFPACK_ZNAME (free_numeric) (&Numeric);
05594 #else
05595 (*current_liboctave_error_handler) ("UMFPACK not installed");
05596 #endif
05597
05598 return Numeric;
05599 }
05600
05601 ComplexMatrix
05602 SparseComplexMatrix::fsolve (MatrixType &mattype, const Matrix& b,
05603 octave_idx_type& err, double& rcond,
05604 solve_singularity_handler sing_handler,
05605 bool calc_cond) const
05606 {
05607 ComplexMatrix retval;
05608
05609 octave_idx_type nr = rows ();
05610 octave_idx_type nc = cols ();
05611 err = 0;
05612
05613 if (nr != nc || nr != b.rows ())
05614 (*current_liboctave_error_handler)
05615 ("matrix dimension mismatch solution of linear equations");
05616 else if (nr == 0 || b.cols () == 0)
05617 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
05618 else
05619 {
05620
05621 volatile int typ = mattype.type ();
05622 mattype.info ();
05623
05624 if (typ == MatrixType::Hermitian)
05625 {
05626 #ifdef HAVE_CHOLMOD
05627 cholmod_common Common;
05628 cholmod_common *cm = &Common;
05629
05630
05631 CHOLMOD_NAME(start) (cm);
05632 cm->prefer_zomplex = false;
05633
05634 double spu = octave_sparse_params::get_key ("spumoni");
05635 if (spu == 0.)
05636 {
05637 cm->print = -1;
05638 cm->print_function = 0;
05639 }
05640 else
05641 {
05642 cm->print = static_cast<int> (spu) + 2;
05643 cm->print_function =&SparseCholPrint;
05644 }
05645
05646 cm->error_handler = &SparseCholError;
05647 cm->complex_divide = CHOLMOD_NAME(divcomplex);
05648 cm->hypotenuse = CHOLMOD_NAME(hypot);
05649
05650 cm->final_ll = true;
05651
05652 cholmod_sparse Astore;
05653 cholmod_sparse *A = &Astore;
05654 double dummy;
05655 A->nrow = nr;
05656 A->ncol = nc;
05657
05658 A->p = cidx();
05659 A->i = ridx();
05660 A->nzmax = nnz();
05661 A->packed = true;
05662 A->sorted = true;
05663 A->nz = 0;
05664 #ifdef IDX_TYPE_LONG
05665 A->itype = CHOLMOD_LONG;
05666 #else
05667 A->itype = CHOLMOD_INT;
05668 #endif
05669 A->dtype = CHOLMOD_DOUBLE;
05670 A->stype = 1;
05671 A->xtype = CHOLMOD_COMPLEX;
05672
05673 if (nr < 1)
05674 A->x = &dummy;
05675 else
05676 A->x = data();
05677
05678 cholmod_dense Bstore;
05679 cholmod_dense *B = &Bstore;
05680 B->nrow = b.rows();
05681 B->ncol = b.cols();
05682 B->d = B->nrow;
05683 B->nzmax = B->nrow * B->ncol;
05684 B->dtype = CHOLMOD_DOUBLE;
05685 B->xtype = CHOLMOD_REAL;
05686 if (nc < 1 || b.cols() < 1)
05687 B->x = &dummy;
05688 else
05689
05690 B->x = const_cast<double *>(b.fortran_vec());
05691
05692 cholmod_factor *L;
05693 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05694 L = CHOLMOD_NAME(analyze) (A, cm);
05695 CHOLMOD_NAME(factorize) (A, L, cm);
05696 if (calc_cond)
05697 rcond = CHOLMOD_NAME(rcond)(L, cm);
05698 else
05699 rcond = 1.;
05700 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05701
05702 if (rcond == 0.0)
05703 {
05704
05705 mattype.mark_as_unsymmetric ();
05706 typ = MatrixType::Full;
05707 }
05708 else
05709 {
05710 volatile double rcond_plus_one = rcond + 1.0;
05711
05712 if (rcond_plus_one == 1.0 || xisnan (rcond))
05713 {
05714 err = -2;
05715
05716 if (sing_handler)
05717 {
05718 sing_handler (rcond);
05719 mattype.mark_as_rectangular ();
05720 }
05721 else
05722 (*current_liboctave_error_handler)
05723 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
05724 rcond);
05725
05726 return retval;
05727 }
05728
05729 cholmod_dense *X;
05730 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05731 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
05732 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05733
05734 retval.resize (b.rows (), b.cols());
05735 for (octave_idx_type j = 0; j < b.cols(); j++)
05736 {
05737 octave_idx_type jr = j * b.rows();
05738 for (octave_idx_type i = 0; i < b.rows(); i++)
05739 retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i];
05740 }
05741
05742 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05743 CHOLMOD_NAME(free_dense) (&X, cm);
05744 CHOLMOD_NAME(free_factor) (&L, cm);
05745 CHOLMOD_NAME(finish) (cm);
05746 static char tmp[] = " ";
05747 CHOLMOD_NAME(print_common) (tmp, cm);
05748 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05749 }
05750 #else
05751 (*current_liboctave_warning_handler)
05752 ("CHOLMOD not installed");
05753
05754 mattype.mark_as_unsymmetric ();
05755 typ = MatrixType::Full;
05756 #endif
05757 }
05758
05759 if (typ == MatrixType::Full)
05760 {
05761 #ifdef HAVE_UMFPACK
05762 Matrix Control, Info;
05763 void *Numeric = factorize (err, rcond, Control, Info,
05764 sing_handler, calc_cond);
05765
05766 if (err == 0)
05767 {
05768 octave_idx_type b_nr = b.rows ();
05769 octave_idx_type b_nc = b.cols ();
05770 int status = 0;
05771 double *control = Control.fortran_vec ();
05772 double *info = Info.fortran_vec ();
05773 const octave_idx_type *Ap = cidx ();
05774 const octave_idx_type *Ai = ridx ();
05775 const Complex *Ax = data ();
05776 #ifdef UMFPACK_SEPARATE_SPLIT
05777 const double *Bx = b.fortran_vec ();
05778 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
05779 for (octave_idx_type i = 0; i < b_nr; i++)
05780 Bz[i] = 0.;
05781 #else
05782 OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
05783 #endif
05784 retval.resize (b_nr, b_nc);
05785 Complex *Xx = retval.fortran_vec ();
05786
05787 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
05788 {
05789 #ifdef UMFPACK_SEPARATE_SPLIT
05790 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
05791 Ai,
05792 reinterpret_cast<const double *> (Ax),
05793 0,
05794 reinterpret_cast<double *> (&Xx[iidx]),
05795 0,
05796 &Bx[iidx], Bz, Numeric,
05797 control, info);
05798 #else
05799 for (octave_idx_type i = 0; i < b_nr; i++)
05800 Bz[i] = b.elem (i, j);
05801
05802 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
05803 Ai,
05804 reinterpret_cast<const double *> (Ax),
05805 0,
05806 reinterpret_cast<double *> (&Xx[iidx]),
05807 0,
05808 reinterpret_cast<const double *> (Bz),
05809 0, Numeric,
05810 control, info);
05811 #endif
05812
05813 if (status < 0)
05814 {
05815 (*current_liboctave_error_handler)
05816 ("SparseComplexMatrix::solve solve failed");
05817
05818 UMFPACK_ZNAME (report_status) (control, status);
05819
05820 err = -1;
05821
05822 break;
05823 }
05824 }
05825
05826 UMFPACK_ZNAME (report_info) (control, info);
05827
05828 UMFPACK_ZNAME (free_numeric) (&Numeric);
05829 }
05830 else
05831 mattype.mark_as_rectangular ();
05832
05833 #else
05834 (*current_liboctave_error_handler) ("UMFPACK not installed");
05835 #endif
05836 }
05837 else if (typ != MatrixType::Hermitian)
05838 (*current_liboctave_error_handler) ("incorrect matrix type");
05839 }
05840
05841 return retval;
05842 }
05843
05844 SparseComplexMatrix
05845 SparseComplexMatrix::fsolve (MatrixType &mattype, const SparseMatrix& b,
05846 octave_idx_type& err, double& rcond,
05847 solve_singularity_handler sing_handler,
05848 bool calc_cond) const
05849 {
05850 SparseComplexMatrix retval;
05851
05852 octave_idx_type nr = rows ();
05853 octave_idx_type nc = cols ();
05854 err = 0;
05855
05856 if (nr != nc || nr != b.rows ())
05857 (*current_liboctave_error_handler)
05858 ("matrix dimension mismatch solution of linear equations");
05859 else if (nr == 0 || b.cols () == 0)
05860 retval = SparseComplexMatrix (nc, b.cols ());
05861 else
05862 {
05863
05864 volatile int typ = mattype.type ();
05865 mattype.info ();
05866
05867 if (typ == MatrixType::Hermitian)
05868 {
05869 #ifdef HAVE_CHOLMOD
05870 cholmod_common Common;
05871 cholmod_common *cm = &Common;
05872
05873
05874 CHOLMOD_NAME(start) (cm);
05875 cm->prefer_zomplex = false;
05876
05877 double spu = octave_sparse_params::get_key ("spumoni");
05878 if (spu == 0.)
05879 {
05880 cm->print = -1;
05881 cm->print_function = 0;
05882 }
05883 else
05884 {
05885 cm->print = static_cast<int> (spu) + 2;
05886 cm->print_function =&SparseCholPrint;
05887 }
05888
05889 cm->error_handler = &SparseCholError;
05890 cm->complex_divide = CHOLMOD_NAME(divcomplex);
05891 cm->hypotenuse = CHOLMOD_NAME(hypot);
05892
05893 cm->final_ll = true;
05894
05895 cholmod_sparse Astore;
05896 cholmod_sparse *A = &Astore;
05897 double dummy;
05898 A->nrow = nr;
05899 A->ncol = nc;
05900
05901 A->p = cidx();
05902 A->i = ridx();
05903 A->nzmax = nnz();
05904 A->packed = true;
05905 A->sorted = true;
05906 A->nz = 0;
05907 #ifdef IDX_TYPE_LONG
05908 A->itype = CHOLMOD_LONG;
05909 #else
05910 A->itype = CHOLMOD_INT;
05911 #endif
05912 A->dtype = CHOLMOD_DOUBLE;
05913 A->stype = 1;
05914 A->xtype = CHOLMOD_COMPLEX;
05915
05916 if (nr < 1)
05917 A->x = &dummy;
05918 else
05919 A->x = data();
05920
05921 cholmod_sparse Bstore;
05922 cholmod_sparse *B = &Bstore;
05923 B->nrow = b.rows();
05924 B->ncol = b.cols();
05925 B->p = b.cidx();
05926 B->i = b.ridx();
05927 B->nzmax = b.nnz();
05928 B->packed = true;
05929 B->sorted = true;
05930 B->nz = 0;
05931 #ifdef IDX_TYPE_LONG
05932 B->itype = CHOLMOD_LONG;
05933 #else
05934 B->itype = CHOLMOD_INT;
05935 #endif
05936 B->dtype = CHOLMOD_DOUBLE;
05937 B->stype = 0;
05938 B->xtype = CHOLMOD_REAL;
05939
05940 if (b.rows() < 1 || b.cols() < 1)
05941 B->x = &dummy;
05942 else
05943 B->x = b.data();
05944
05945 cholmod_factor *L;
05946 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05947 L = CHOLMOD_NAME(analyze) (A, cm);
05948 CHOLMOD_NAME(factorize) (A, L, cm);
05949 if (calc_cond)
05950 rcond = CHOLMOD_NAME(rcond)(L, cm);
05951 else
05952 rcond = 1.;
05953 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05954
05955 if (rcond == 0.0)
05956 {
05957
05958 mattype.mark_as_unsymmetric ();
05959 typ = MatrixType::Full;
05960 }
05961 else
05962 {
05963 volatile double rcond_plus_one = rcond + 1.0;
05964
05965 if (rcond_plus_one == 1.0 || xisnan (rcond))
05966 {
05967 err = -2;
05968
05969 if (sing_handler)
05970 {
05971 sing_handler (rcond);
05972 mattype.mark_as_rectangular ();
05973 }
05974 else
05975 (*current_liboctave_error_handler)
05976 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
05977 rcond);
05978
05979 return retval;
05980 }
05981
05982 cholmod_sparse *X;
05983 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05984 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
05985 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
05986
05987 retval = SparseComplexMatrix
05988 (static_cast<octave_idx_type>(X->nrow),
05989 static_cast<octave_idx_type>(X->ncol),
05990 static_cast<octave_idx_type>(X->nzmax));
05991 for (octave_idx_type j = 0;
05992 j <= static_cast<octave_idx_type>(X->ncol); j++)
05993 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
05994 for (octave_idx_type j = 0;
05995 j < static_cast<octave_idx_type>(X->nzmax); j++)
05996 {
05997 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
05998 retval.xdata(j) = static_cast<Complex *>(X->x)[j];
05999 }
06000
06001 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06002 CHOLMOD_NAME(free_sparse) (&X, cm);
06003 CHOLMOD_NAME(free_factor) (&L, cm);
06004 CHOLMOD_NAME(finish) (cm);
06005 static char tmp[] = " ";
06006 CHOLMOD_NAME(print_common) (tmp, cm);
06007 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06008 }
06009 #else
06010 (*current_liboctave_warning_handler)
06011 ("CHOLMOD not installed");
06012
06013 mattype.mark_as_unsymmetric ();
06014 typ = MatrixType::Full;
06015 #endif
06016 }
06017
06018 if (typ == MatrixType::Full)
06019 {
06020 #ifdef HAVE_UMFPACK
06021 Matrix Control, Info;
06022 void *Numeric = factorize (err, rcond, Control, Info,
06023 sing_handler, calc_cond);
06024
06025 if (err == 0)
06026 {
06027 octave_idx_type b_nr = b.rows ();
06028 octave_idx_type b_nc = b.cols ();
06029 int status = 0;
06030 double *control = Control.fortran_vec ();
06031 double *info = Info.fortran_vec ();
06032 const octave_idx_type *Ap = cidx ();
06033 const octave_idx_type *Ai = ridx ();
06034 const Complex *Ax = data ();
06035
06036 #ifdef UMFPACK_SEPARATE_SPLIT
06037 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
06038 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
06039 for (octave_idx_type i = 0; i < b_nr; i++)
06040 Bz[i] = 0.;
06041 #else
06042 OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
06043 #endif
06044
06045
06046
06047 octave_idx_type x_nz = b.nnz ();
06048 octave_idx_type ii = 0;
06049 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
06050
06051 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
06052
06053 retval.xcidx(0) = 0;
06054 for (octave_idx_type j = 0; j < b_nc; j++)
06055 {
06056
06057 #ifdef UMFPACK_SEPARATE_SPLIT
06058 for (octave_idx_type i = 0; i < b_nr; i++)
06059 Bx[i] = b.elem (i, j);
06060
06061 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
06062 Ai,
06063 reinterpret_cast<const double *> (Ax),
06064 0,
06065 reinterpret_cast<double *> (Xx),
06066 0,
06067 Bx, Bz, Numeric, control,
06068 info);
06069 #else
06070 for (octave_idx_type i = 0; i < b_nr; i++)
06071 Bz[i] = b.elem (i, j);
06072
06073 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
06074 reinterpret_cast<const double *> (Ax),
06075 0,
06076 reinterpret_cast<double *> (Xx),
06077 0,
06078 reinterpret_cast<double *> (Bz),
06079 0,
06080 Numeric, control,
06081 info);
06082 #endif
06083 if (status < 0)
06084 {
06085 (*current_liboctave_error_handler)
06086 ("SparseComplexMatrix::solve solve failed");
06087
06088 UMFPACK_ZNAME (report_status) (control, status);
06089
06090 err = -1;
06091
06092 break;
06093 }
06094
06095 for (octave_idx_type i = 0; i < b_nr; i++)
06096 {
06097 Complex tmp = Xx[i];
06098 if (tmp != 0.0)
06099 {
06100 if (ii == x_nz)
06101 {
06102
06103 octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
06104 sz = (sz > 10 ? sz : 10) + x_nz;
06105 retval.change_capacity (sz);
06106 x_nz = sz;
06107 }
06108 retval.xdata(ii) = tmp;
06109 retval.xridx(ii++) = i;
06110 }
06111 }
06112 retval.xcidx(j+1) = ii;
06113 }
06114
06115 retval.maybe_compress ();
06116
06117 UMFPACK_ZNAME (report_info) (control, info);
06118
06119 UMFPACK_ZNAME (free_numeric) (&Numeric);
06120 }
06121 else
06122 mattype.mark_as_rectangular ();
06123
06124 #else
06125 (*current_liboctave_error_handler) ("UMFPACK not installed");
06126 #endif
06127 }
06128 else if (typ != MatrixType::Hermitian)
06129 (*current_liboctave_error_handler) ("incorrect matrix type");
06130 }
06131
06132 return retval;
06133 }
06134
06135 ComplexMatrix
06136 SparseComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b,
06137 octave_idx_type& err, double& rcond,
06138 solve_singularity_handler sing_handler,
06139 bool calc_cond) const
06140 {
06141 ComplexMatrix retval;
06142
06143 octave_idx_type nr = rows ();
06144 octave_idx_type nc = cols ();
06145 err = 0;
06146
06147 if (nr != nc || nr != b.rows ())
06148 (*current_liboctave_error_handler)
06149 ("matrix dimension mismatch solution of linear equations");
06150 else if (nr == 0 || b.cols () == 0)
06151 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
06152 else
06153 {
06154
06155 volatile int typ = mattype.type ();
06156 mattype.info ();
06157
06158 if (typ == MatrixType::Hermitian)
06159 {
06160 #ifdef HAVE_CHOLMOD
06161 cholmod_common Common;
06162 cholmod_common *cm = &Common;
06163
06164
06165 CHOLMOD_NAME(start) (cm);
06166 cm->prefer_zomplex = false;
06167
06168 double spu = octave_sparse_params::get_key ("spumoni");
06169 if (spu == 0.)
06170 {
06171 cm->print = -1;
06172 cm->print_function = 0;
06173 }
06174 else
06175 {
06176 cm->print = static_cast<int> (spu) + 2;
06177 cm->print_function =&SparseCholPrint;
06178 }
06179
06180 cm->error_handler = &SparseCholError;
06181 cm->complex_divide = CHOLMOD_NAME(divcomplex);
06182 cm->hypotenuse = CHOLMOD_NAME(hypot);
06183
06184 cm->final_ll = true;
06185
06186 cholmod_sparse Astore;
06187 cholmod_sparse *A = &Astore;
06188 double dummy;
06189 A->nrow = nr;
06190 A->ncol = nc;
06191
06192 A->p = cidx();
06193 A->i = ridx();
06194 A->nzmax = nnz();
06195 A->packed = true;
06196 A->sorted = true;
06197 A->nz = 0;
06198 #ifdef IDX_TYPE_LONG
06199 A->itype = CHOLMOD_LONG;
06200 #else
06201 A->itype = CHOLMOD_INT;
06202 #endif
06203 A->dtype = CHOLMOD_DOUBLE;
06204 A->stype = 1;
06205 A->xtype = CHOLMOD_COMPLEX;
06206
06207 if (nr < 1)
06208 A->x = &dummy;
06209 else
06210 A->x = data();
06211
06212 cholmod_dense Bstore;
06213 cholmod_dense *B = &Bstore;
06214 B->nrow = b.rows();
06215 B->ncol = b.cols();
06216 B->d = B->nrow;
06217 B->nzmax = B->nrow * B->ncol;
06218 B->dtype = CHOLMOD_DOUBLE;
06219 B->xtype = CHOLMOD_COMPLEX;
06220 if (nc < 1 || b.cols() < 1)
06221 B->x = &dummy;
06222 else
06223
06224 B->x = const_cast<Complex *>(b.fortran_vec());
06225
06226 cholmod_factor *L;
06227 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06228 L = CHOLMOD_NAME(analyze) (A, cm);
06229 CHOLMOD_NAME(factorize) (A, L, cm);
06230 if (calc_cond)
06231 rcond = CHOLMOD_NAME(rcond)(L, cm);
06232 else
06233 rcond = 1.;
06234 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06235
06236 if (rcond == 0.0)
06237 {
06238
06239 mattype.mark_as_unsymmetric ();
06240 typ = MatrixType::Full;
06241 }
06242 else
06243 {
06244 volatile double rcond_plus_one = rcond + 1.0;
06245
06246 if (rcond_plus_one == 1.0 || xisnan (rcond))
06247 {
06248 err = -2;
06249
06250 if (sing_handler)
06251 {
06252 sing_handler (rcond);
06253 mattype.mark_as_rectangular ();
06254 }
06255 else
06256 (*current_liboctave_error_handler)
06257 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
06258 rcond);
06259
06260 return retval;
06261 }
06262
06263 cholmod_dense *X;
06264 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06265 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
06266 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06267
06268 retval.resize (b.rows (), b.cols());
06269 for (octave_idx_type j = 0; j < b.cols(); j++)
06270 {
06271 octave_idx_type jr = j * b.rows();
06272 for (octave_idx_type i = 0; i < b.rows(); i++)
06273 retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i];
06274 }
06275
06276 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06277 CHOLMOD_NAME(free_dense) (&X, cm);
06278 CHOLMOD_NAME(free_factor) (&L, cm);
06279 CHOLMOD_NAME(finish) (cm);
06280 static char tmp[] = " ";
06281 CHOLMOD_NAME(print_common) (tmp, cm);
06282 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06283 }
06284 #else
06285 (*current_liboctave_warning_handler)
06286 ("CHOLMOD not installed");
06287
06288 mattype.mark_as_unsymmetric ();
06289 typ = MatrixType::Full;
06290 #endif
06291 }
06292
06293 if (typ == MatrixType::Full)
06294 {
06295 #ifdef HAVE_UMFPACK
06296 Matrix Control, Info;
06297 void *Numeric = factorize (err, rcond, Control, Info,
06298 sing_handler, calc_cond);
06299
06300 if (err == 0)
06301 {
06302 octave_idx_type b_nr = b.rows ();
06303 octave_idx_type b_nc = b.cols ();
06304 int status = 0;
06305 double *control = Control.fortran_vec ();
06306 double *info = Info.fortran_vec ();
06307 const octave_idx_type *Ap = cidx ();
06308 const octave_idx_type *Ai = ridx ();
06309 const Complex *Ax = data ();
06310 const Complex *Bx = b.fortran_vec ();
06311
06312 retval.resize (b_nr, b_nc);
06313 Complex *Xx = retval.fortran_vec ();
06314
06315 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
06316 {
06317 status =
06318 UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
06319 reinterpret_cast<const double *> (Ax),
06320 0,
06321 reinterpret_cast<double *> (&Xx[iidx]),
06322 0,
06323 reinterpret_cast<const double *> (&Bx[iidx]),
06324 0, Numeric, control, info);
06325
06326 if (status < 0)
06327 {
06328 (*current_liboctave_error_handler)
06329 ("SparseComplexMatrix::solve solve failed");
06330
06331 UMFPACK_ZNAME (report_status) (control, status);
06332
06333 err = -1;
06334
06335 break;
06336 }
06337 }
06338
06339 UMFPACK_ZNAME (report_info) (control, info);
06340
06341 UMFPACK_ZNAME (free_numeric) (&Numeric);
06342 }
06343 else
06344 mattype.mark_as_rectangular ();
06345
06346 #else
06347 (*current_liboctave_error_handler) ("UMFPACK not installed");
06348 #endif
06349 }
06350 else if (typ != MatrixType::Hermitian)
06351 (*current_liboctave_error_handler) ("incorrect matrix type");
06352 }
06353
06354 return retval;
06355 }
06356
06357 SparseComplexMatrix
06358 SparseComplexMatrix::fsolve (MatrixType &mattype, const SparseComplexMatrix& b,
06359 octave_idx_type& err, double& rcond,
06360 solve_singularity_handler sing_handler,
06361 bool calc_cond) const
06362 {
06363 SparseComplexMatrix retval;
06364
06365 octave_idx_type nr = rows ();
06366 octave_idx_type nc = cols ();
06367 err = 0;
06368
06369 if (nr != nc || nr != b.rows ())
06370 (*current_liboctave_error_handler)
06371 ("matrix dimension mismatch solution of linear equations");
06372 else if (nr == 0 || b.cols () == 0)
06373 retval = SparseComplexMatrix (nc, b.cols ());
06374 else
06375 {
06376
06377 volatile int typ = mattype.type ();
06378 mattype.info ();
06379
06380 if (typ == MatrixType::Hermitian)
06381 {
06382 #ifdef HAVE_CHOLMOD
06383 cholmod_common Common;
06384 cholmod_common *cm = &Common;
06385
06386
06387 CHOLMOD_NAME(start) (cm);
06388 cm->prefer_zomplex = false;
06389
06390 double spu = octave_sparse_params::get_key ("spumoni");
06391 if (spu == 0.)
06392 {
06393 cm->print = -1;
06394 cm->print_function = 0;
06395 }
06396 else
06397 {
06398 cm->print = static_cast<int> (spu) + 2;
06399 cm->print_function =&SparseCholPrint;
06400 }
06401
06402 cm->error_handler = &SparseCholError;
06403 cm->complex_divide = CHOLMOD_NAME(divcomplex);
06404 cm->hypotenuse = CHOLMOD_NAME(hypot);
06405
06406 cm->final_ll = true;
06407
06408 cholmod_sparse Astore;
06409 cholmod_sparse *A = &Astore;
06410 double dummy;
06411 A->nrow = nr;
06412 A->ncol = nc;
06413
06414 A->p = cidx();
06415 A->i = ridx();
06416 A->nzmax = nnz();
06417 A->packed = true;
06418 A->sorted = true;
06419 A->nz = 0;
06420 #ifdef IDX_TYPE_LONG
06421 A->itype = CHOLMOD_LONG;
06422 #else
06423 A->itype = CHOLMOD_INT;
06424 #endif
06425 A->dtype = CHOLMOD_DOUBLE;
06426 A->stype = 1;
06427 A->xtype = CHOLMOD_COMPLEX;
06428
06429 if (nr < 1)
06430 A->x = &dummy;
06431 else
06432 A->x = data();
06433
06434 cholmod_sparse Bstore;
06435 cholmod_sparse *B = &Bstore;
06436 B->nrow = b.rows();
06437 B->ncol = b.cols();
06438 B->p = b.cidx();
06439 B->i = b.ridx();
06440 B->nzmax = b.nnz();
06441 B->packed = true;
06442 B->sorted = true;
06443 B->nz = 0;
06444 #ifdef IDX_TYPE_LONG
06445 B->itype = CHOLMOD_LONG;
06446 #else
06447 B->itype = CHOLMOD_INT;
06448 #endif
06449 B->dtype = CHOLMOD_DOUBLE;
06450 B->stype = 0;
06451 B->xtype = CHOLMOD_COMPLEX;
06452
06453 if (b.rows() < 1 || b.cols() < 1)
06454 B->x = &dummy;
06455 else
06456 B->x = b.data();
06457
06458 cholmod_factor *L;
06459 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06460 L = CHOLMOD_NAME(analyze) (A, cm);
06461 CHOLMOD_NAME(factorize) (A, L, cm);
06462 if (calc_cond)
06463 rcond = CHOLMOD_NAME(rcond)(L, cm);
06464 else
06465 rcond = 1.;
06466 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06467
06468 if (rcond == 0.0)
06469 {
06470
06471 mattype.mark_as_unsymmetric ();
06472 typ = MatrixType::Full;
06473 }
06474 else
06475 {
06476 volatile double rcond_plus_one = rcond + 1.0;
06477
06478 if (rcond_plus_one == 1.0 || xisnan (rcond))
06479 {
06480 err = -2;
06481
06482 if (sing_handler)
06483 {
06484 sing_handler (rcond);
06485 mattype.mark_as_rectangular ();
06486 }
06487 else
06488 (*current_liboctave_error_handler)
06489 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
06490 rcond);
06491
06492 return retval;
06493 }
06494
06495 cholmod_sparse *X;
06496 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06497 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
06498 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06499
06500 retval = SparseComplexMatrix
06501 (static_cast<octave_idx_type>(X->nrow),
06502 static_cast<octave_idx_type>(X->ncol),
06503 static_cast<octave_idx_type>(X->nzmax));
06504 for (octave_idx_type j = 0;
06505 j <= static_cast<octave_idx_type>(X->ncol); j++)
06506 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
06507 for (octave_idx_type j = 0;
06508 j < static_cast<octave_idx_type>(X->nzmax); j++)
06509 {
06510 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
06511 retval.xdata(j) = static_cast<Complex *>(X->x)[j];
06512 }
06513
06514 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06515 CHOLMOD_NAME(free_sparse) (&X, cm);
06516 CHOLMOD_NAME(free_factor) (&L, cm);
06517 CHOLMOD_NAME(finish) (cm);
06518 static char tmp[] = " ";
06519 CHOLMOD_NAME(print_common) (tmp, cm);
06520 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
06521 }
06522 #else
06523 (*current_liboctave_warning_handler)
06524 ("CHOLMOD not installed");
06525
06526 mattype.mark_as_unsymmetric ();
06527 typ = MatrixType::Full;
06528 #endif
06529 }
06530
06531 if (typ == MatrixType::Full)
06532 {
06533 #ifdef HAVE_UMFPACK
06534 Matrix Control, Info;
06535 void *Numeric = factorize (err, rcond, Control, Info,
06536 sing_handler, calc_cond);
06537
06538 if (err == 0)
06539 {
06540 octave_idx_type b_nr = b.rows ();
06541 octave_idx_type b_nc = b.cols ();
06542 int status = 0;
06543 double *control = Control.fortran_vec ();
06544 double *info = Info.fortran_vec ();
06545 const octave_idx_type *Ap = cidx ();
06546 const octave_idx_type *Ai = ridx ();
06547 const Complex *Ax = data ();
06548
06549 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
06550
06551
06552
06553 octave_idx_type x_nz = b.nnz ();
06554 octave_idx_type ii = 0;
06555 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
06556
06557 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
06558
06559 retval.xcidx(0) = 0;
06560 for (octave_idx_type j = 0; j < b_nc; j++)
06561 {
06562 for (octave_idx_type i = 0; i < b_nr; i++)
06563 Bx[i] = b (i,j);
06564
06565 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
06566 Ai,
06567 reinterpret_cast<const double *> (Ax),
06568 0,
06569 reinterpret_cast<double *> (Xx),
06570 0,
06571 reinterpret_cast<double *> (Bx),
06572 0, Numeric, control, info);
06573
06574 if (status < 0)
06575 {
06576 (*current_liboctave_error_handler)
06577 ("SparseComplexMatrix::solve solve failed");
06578
06579 UMFPACK_ZNAME (report_status) (control, status);
06580
06581 err = -1;
06582
06583 break;
06584 }
06585
06586 for (octave_idx_type i = 0; i < b_nr; i++)
06587 {
06588 Complex tmp = Xx[i];
06589 if (tmp != 0.0)
06590 {
06591 if (ii == x_nz)
06592 {
06593
06594 octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
06595 sz = (sz > 10 ? sz : 10) + x_nz;
06596 retval.change_capacity (sz);
06597 x_nz = sz;
06598 }
06599 retval.xdata(ii) = tmp;
06600 retval.xridx(ii++) = i;
06601 }
06602 }
06603 retval.xcidx(j+1) = ii;
06604 }
06605
06606 retval.maybe_compress ();
06607
06608 rcond = Info (UMFPACK_RCOND);
06609 volatile double rcond_plus_one = rcond + 1.0;
06610
06611 if (status == UMFPACK_WARNING_singular_matrix ||
06612 rcond_plus_one == 1.0 || xisnan (rcond))
06613 {
06614 err = -2;
06615
06616 if (sing_handler)
06617 sing_handler (rcond);
06618 else
06619 (*current_liboctave_error_handler)
06620 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
06621 rcond);
06622
06623 }
06624
06625 UMFPACK_ZNAME (report_info) (control, info);
06626
06627 UMFPACK_ZNAME (free_numeric) (&Numeric);
06628 }
06629 else
06630 mattype.mark_as_rectangular ();
06631
06632 #else
06633 (*current_liboctave_error_handler) ("UMFPACK not installed");
06634 #endif
06635 }
06636 else if (typ != MatrixType::Hermitian)
06637 (*current_liboctave_error_handler) ("incorrect matrix type");
06638 }
06639
06640 return retval;
06641 }
06642
06643 ComplexMatrix
06644 SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b) const
06645 {
06646 octave_idx_type info;
06647 double rcond;
06648 return solve (mattype, b, info, rcond, 0);
06649 }
06650
06651 ComplexMatrix
06652 SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b,
06653 octave_idx_type& info) const
06654 {
06655 double rcond;
06656 return solve (mattype, b, info, rcond, 0);
06657 }
06658
06659 ComplexMatrix
06660 SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b,
06661 octave_idx_type& info, double& rcond) const
06662 {
06663 return solve (mattype, b, info, rcond, 0);
06664 }
06665
06666 ComplexMatrix
06667 SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b,
06668 octave_idx_type& err, double& rcond,
06669 solve_singularity_handler sing_handler,
06670 bool singular_fallback) const
06671 {
06672 ComplexMatrix retval;
06673 int typ = mattype.type (false);
06674
06675 if (typ == MatrixType::Unknown)
06676 typ = mattype.type (*this);
06677
06678 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
06679 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
06680 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
06681 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
06682 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
06683 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
06684 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
06685 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
06686 else if (typ == MatrixType::Tridiagonal ||
06687 typ == MatrixType::Tridiagonal_Hermitian)
06688 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
06689 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
06690 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
06691 else if (typ != MatrixType::Rectangular)
06692 {
06693 (*current_liboctave_error_handler) ("unknown matrix type");
06694 return ComplexMatrix ();
06695 }
06696
06697 if (singular_fallback && mattype.type(false) == MatrixType::Rectangular)
06698 {
06699 rcond = 1.;
06700 #ifdef USE_QRSOLVE
06701 retval = qrsolve (*this, b, err);
06702 #else
06703 retval = dmsolve<ComplexMatrix, SparseComplexMatrix,
06704 Matrix> (*this, b, err);
06705 #endif
06706 }
06707
06708 return retval;
06709 }
06710
06711 SparseComplexMatrix
06712 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b) const
06713 {
06714 octave_idx_type info;
06715 double rcond;
06716 return solve (mattype, b, info, rcond, 0);
06717 }
06718
06719 SparseComplexMatrix
06720 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
06721 octave_idx_type& info) const
06722 {
06723 double rcond;
06724 return solve (mattype, b, info, rcond, 0);
06725 }
06726
06727 SparseComplexMatrix
06728 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
06729 octave_idx_type& info, double& rcond) const
06730 {
06731 return solve (mattype, b, info, rcond, 0);
06732 }
06733
06734 SparseComplexMatrix
06735 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
06736 octave_idx_type& err, double& rcond,
06737 solve_singularity_handler sing_handler,
06738 bool singular_fallback) const
06739 {
06740 SparseComplexMatrix retval;
06741 int typ = mattype.type (false);
06742
06743 if (typ == MatrixType::Unknown)
06744 typ = mattype.type (*this);
06745
06746 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
06747 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
06748 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
06749 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
06750 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
06751 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
06752 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
06753 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
06754 else if (typ == MatrixType::Tridiagonal ||
06755 typ == MatrixType::Tridiagonal_Hermitian)
06756 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
06757 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
06758 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
06759 else if (typ != MatrixType::Rectangular)
06760 {
06761 (*current_liboctave_error_handler) ("unknown matrix type");
06762 return SparseComplexMatrix ();
06763 }
06764
06765 if (singular_fallback && mattype.type(false) == MatrixType::Rectangular)
06766 {
06767 rcond = 1.;
06768 #ifdef USE_QRSOLVE
06769 retval = qrsolve (*this, b, err);
06770 #else
06771 retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix,
06772 SparseMatrix> (*this, b, err);
06773 #endif
06774 }
06775
06776 return retval;
06777 }
06778
06779 ComplexMatrix
06780 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b) const
06781 {
06782 octave_idx_type info;
06783 double rcond;
06784 return solve (mattype, b, info, rcond, 0);
06785 }
06786
06787 ComplexMatrix
06788 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
06789 octave_idx_type& info) const
06790 {
06791 double rcond;
06792 return solve (mattype, b, info, rcond, 0);
06793 }
06794
06795 ComplexMatrix
06796 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
06797 octave_idx_type& info, double& rcond) const
06798 {
06799 return solve (mattype, b, info, rcond, 0);
06800 }
06801
06802 ComplexMatrix
06803 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
06804 octave_idx_type& err, double& rcond,
06805 solve_singularity_handler sing_handler,
06806 bool singular_fallback) const
06807 {
06808 ComplexMatrix retval;
06809 int typ = mattype.type (false);
06810
06811 if (typ == MatrixType::Unknown)
06812 typ = mattype.type (*this);
06813
06814 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
06815 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
06816 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
06817 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
06818 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
06819 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
06820 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
06821 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
06822 else if (typ == MatrixType::Tridiagonal ||
06823 typ == MatrixType::Tridiagonal_Hermitian)
06824 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
06825 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
06826 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
06827 else if (typ != MatrixType::Rectangular)
06828 {
06829 (*current_liboctave_error_handler) ("unknown matrix type");
06830 return ComplexMatrix ();
06831 }
06832
06833 if (singular_fallback && mattype.type(false) == MatrixType::Rectangular)
06834 {
06835 rcond = 1.;
06836 #ifdef USE_QRSOLVE
06837 retval = qrsolve (*this, b, err);
06838 #else
06839 retval = dmsolve<ComplexMatrix, SparseComplexMatrix,
06840 ComplexMatrix> (*this, b, err);
06841 #endif
06842 }
06843
06844 return retval;
06845 }
06846
06847 SparseComplexMatrix
06848 SparseComplexMatrix::solve (MatrixType &mattype,
06849 const SparseComplexMatrix& b) const
06850 {
06851 octave_idx_type info;
06852 double rcond;
06853 return solve (mattype, b, info, rcond, 0);
06854 }
06855
06856 SparseComplexMatrix
06857 SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
06858 octave_idx_type& info) const
06859 {
06860 double rcond;
06861 return solve (mattype, b, info, rcond, 0);
06862 }
06863
06864 SparseComplexMatrix
06865 SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
06866 octave_idx_type& info, double& rcond) const
06867 {
06868 return solve (mattype, b, info, rcond, 0);
06869 }
06870
06871 SparseComplexMatrix
06872 SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
06873 octave_idx_type& err, double& rcond,
06874 solve_singularity_handler sing_handler,
06875 bool singular_fallback) const
06876 {
06877 SparseComplexMatrix retval;
06878 int typ = mattype.type (false);
06879
06880 if (typ == MatrixType::Unknown)
06881 typ = mattype.type (*this);
06882
06883 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
06884 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
06885 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
06886 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
06887 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
06888 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
06889 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
06890 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
06891 else if (typ == MatrixType::Tridiagonal ||
06892 typ == MatrixType::Tridiagonal_Hermitian)
06893 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
06894 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
06895 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
06896 else if (typ != MatrixType::Rectangular)
06897 {
06898 (*current_liboctave_error_handler) ("unknown matrix type");
06899 return SparseComplexMatrix ();
06900 }
06901
06902 if (singular_fallback && mattype.type(false) == MatrixType::Rectangular)
06903 {
06904 rcond = 1.;
06905 #ifdef USE_QRSOLVE
06906 retval = qrsolve (*this, b, err);
06907 #else
06908 retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix,
06909 SparseComplexMatrix> (*this, b, err);
06910 #endif
06911 }
06912
06913 return retval;
06914 }
06915
06916 ComplexColumnVector
06917 SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b) const
06918 {
06919 octave_idx_type info; double rcond;
06920 return solve (mattype, b, info, rcond);
06921 }
06922
06923 ComplexColumnVector
06924 SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b,
06925 octave_idx_type& info) const
06926 {
06927 double rcond;
06928 return solve (mattype, b, info, rcond);
06929 }
06930
06931 ComplexColumnVector
06932 SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b,
06933 octave_idx_type& info, double& rcond) const
06934 {
06935 return solve (mattype, b, info, rcond, 0);
06936 }
06937
06938 ComplexColumnVector
06939 SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b,
06940 octave_idx_type& info, double& rcond,
06941 solve_singularity_handler sing_handler) const
06942 {
06943 Matrix tmp (b);
06944 return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
06945 }
06946
06947 ComplexColumnVector
06948 SparseComplexMatrix::solve (MatrixType &mattype,
06949 const ComplexColumnVector& b) const
06950 {
06951 octave_idx_type info;
06952 double rcond;
06953 return solve (mattype, b, info, rcond, 0);
06954 }
06955
06956 ComplexColumnVector
06957 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b,
06958 octave_idx_type& info) const
06959 {
06960 double rcond;
06961 return solve (mattype, b, info, rcond, 0);
06962 }
06963
06964 ComplexColumnVector
06965 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b,
06966 octave_idx_type& info, double& rcond) const
06967 {
06968 return solve (mattype, b, info, rcond, 0);
06969 }
06970
06971 ComplexColumnVector
06972 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b,
06973 octave_idx_type& info, double& rcond,
06974 solve_singularity_handler sing_handler) const
06975 {
06976 ComplexMatrix tmp (b);
06977 return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
06978 }
06979
06980 ComplexMatrix
06981 SparseComplexMatrix::solve (const Matrix& b) const
06982 {
06983 octave_idx_type info;
06984 double rcond;
06985 return solve (b, info, rcond, 0);
06986 }
06987
06988 ComplexMatrix
06989 SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const
06990 {
06991 double rcond;
06992 return solve (b, info, rcond, 0);
06993 }
06994
06995 ComplexMatrix
06996 SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info,
06997 double& rcond) const
06998 {
06999 return solve (b, info, rcond, 0);
07000 }
07001
07002 ComplexMatrix
07003 SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& err,
07004 double& rcond,
07005 solve_singularity_handler sing_handler) const
07006 {
07007 MatrixType mattype (*this);
07008 return solve (mattype, b, err, rcond, sing_handler);
07009 }
07010
07011 SparseComplexMatrix
07012 SparseComplexMatrix::solve (const SparseMatrix& b) const
07013 {
07014 octave_idx_type info;
07015 double rcond;
07016 return solve (b, info, rcond, 0);
07017 }
07018
07019 SparseComplexMatrix
07020 SparseComplexMatrix::solve (const SparseMatrix& b,
07021 octave_idx_type& info) const
07022 {
07023 double rcond;
07024 return solve (b, info, rcond, 0);
07025 }
07026
07027 SparseComplexMatrix
07028 SparseComplexMatrix::solve (const SparseMatrix& b,
07029 octave_idx_type& info, double& rcond) const
07030 {
07031 return solve (b, info, rcond, 0);
07032 }
07033
07034 SparseComplexMatrix
07035 SparseComplexMatrix::solve (const SparseMatrix& b,
07036 octave_idx_type& err, double& rcond,
07037 solve_singularity_handler sing_handler) const
07038 {
07039 MatrixType mattype (*this);
07040 return solve (mattype, b, err, rcond, sing_handler);
07041 }
07042
07043 ComplexMatrix
07044 SparseComplexMatrix::solve (const ComplexMatrix& b,
07045 octave_idx_type& info) const
07046 {
07047 double rcond;
07048 return solve (b, info, rcond, 0);
07049 }
07050
07051 ComplexMatrix
07052 SparseComplexMatrix::solve (const ComplexMatrix& b,
07053 octave_idx_type& info, double& rcond) const
07054 {
07055 return solve (b, info, rcond, 0);
07056 }
07057
07058 ComplexMatrix
07059 SparseComplexMatrix::solve (const ComplexMatrix& b,
07060 octave_idx_type& err, double& rcond,
07061 solve_singularity_handler sing_handler) const
07062 {
07063 MatrixType mattype (*this);
07064 return solve (mattype, b, err, rcond, sing_handler);
07065 }
07066
07067 SparseComplexMatrix
07068 SparseComplexMatrix::solve (const SparseComplexMatrix& b) const
07069 {
07070 octave_idx_type info;
07071 double rcond;
07072 return solve (b, info, rcond, 0);
07073 }
07074
07075 SparseComplexMatrix
07076 SparseComplexMatrix::solve (const SparseComplexMatrix& b,
07077 octave_idx_type& info) const
07078 {
07079 double rcond;
07080 return solve (b, info, rcond, 0);
07081 }
07082
07083 SparseComplexMatrix
07084 SparseComplexMatrix::solve (const SparseComplexMatrix& b,
07085 octave_idx_type& info, double& rcond) const
07086 {
07087 return solve (b, info, rcond, 0);
07088 }
07089
07090 SparseComplexMatrix
07091 SparseComplexMatrix::solve (const SparseComplexMatrix& b,
07092 octave_idx_type& err, double& rcond,
07093 solve_singularity_handler sing_handler) const
07094 {
07095 MatrixType mattype (*this);
07096 return solve (mattype, b, err, rcond, sing_handler);
07097 }
07098
07099 ComplexColumnVector
07100 SparseComplexMatrix::solve (const ColumnVector& b) const
07101 {
07102 octave_idx_type info; double rcond;
07103 return solve (b, info, rcond);
07104 }
07105
07106 ComplexColumnVector
07107 SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info) const
07108 {
07109 double rcond;
07110 return solve (b, info, rcond);
07111 }
07112
07113 ComplexColumnVector
07114 SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
07115 double& rcond) const
07116 {
07117 return solve (b, info, rcond, 0);
07118 }
07119
07120 ComplexColumnVector
07121 SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond,
07122 solve_singularity_handler sing_handler) const
07123 {
07124 Matrix tmp (b);
07125 return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
07126 }
07127
07128 ComplexColumnVector
07129 SparseComplexMatrix::solve (const ComplexColumnVector& b) const
07130 {
07131 octave_idx_type info;
07132 double rcond;
07133 return solve (b, info, rcond, 0);
07134 }
07135
07136 ComplexColumnVector
07137 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const
07138 {
07139 double rcond;
07140 return solve (b, info, rcond, 0);
07141 }
07142
07143 ComplexColumnVector
07144 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
07145 double& rcond) const
07146 {
07147 return solve (b, info, rcond, 0);
07148 }
07149
07150 ComplexColumnVector
07151 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
07152 double& rcond,
07153 solve_singularity_handler sing_handler) const
07154 {
07155 ComplexMatrix tmp (b);
07156 return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
07157 }
07158
07159
07160 SparseBoolMatrix
07161 SparseComplexMatrix::operator ! (void) const
07162 {
07163 if (any_element_is_nan ())
07164 gripe_nan_to_logical_conversion ();
07165
07166 octave_idx_type nr = rows ();
07167 octave_idx_type nc = cols ();
07168 octave_idx_type nz1 = nnz ();
07169 octave_idx_type nz2 = nr*nc - nz1;
07170
07171 SparseBoolMatrix r (nr, nc, nz2);
07172
07173 octave_idx_type ii = 0;
07174 octave_idx_type jj = 0;
07175 r.cidx (0) = 0;
07176 for (octave_idx_type i = 0; i < nc; i++)
07177 {
07178 for (octave_idx_type j = 0; j < nr; j++)
07179 {
07180 if (jj < cidx(i+1) && ridx(jj) == j)
07181 jj++;
07182 else
07183 {
07184 r.data(ii) = true;
07185 r.ridx(ii++) = j;
07186 }
07187 }
07188 r.cidx (i+1) = ii;
07189 }
07190
07191 return r;
07192 }
07193
07194 SparseComplexMatrix
07195 SparseComplexMatrix::squeeze (void) const
07196 {
07197 return MSparse<Complex>::squeeze ();
07198 }
07199
07200 SparseComplexMatrix
07201 SparseComplexMatrix::reshape (const dim_vector& new_dims) const
07202 {
07203 return MSparse<Complex>::reshape (new_dims);
07204 }
07205
07206 SparseComplexMatrix
07207 SparseComplexMatrix::permute (const Array<octave_idx_type>& vec, bool inv) const
07208 {
07209 return MSparse<Complex>::permute (vec, inv);
07210 }
07211
07212 SparseComplexMatrix
07213 SparseComplexMatrix::ipermute (const Array<octave_idx_type>& vec) const
07214 {
07215 return MSparse<Complex>::ipermute (vec);
07216 }
07217
07218
07219
07220 bool
07221 SparseComplexMatrix::any_element_is_nan (void) const
07222 {
07223 octave_idx_type nel = nnz ();
07224
07225 for (octave_idx_type i = 0; i < nel; i++)
07226 {
07227 Complex val = data (i);
07228 if (xisnan (val))
07229 return true;
07230 }
07231
07232 return false;
07233 }
07234
07235 bool
07236 SparseComplexMatrix::any_element_is_inf_or_nan (void) const
07237 {
07238 octave_idx_type nel = nnz ();
07239
07240 for (octave_idx_type i = 0; i < nel; i++)
07241 {
07242 Complex val = data (i);
07243 if (xisinf (val) || xisnan (val))
07244 return true;
07245 }
07246
07247 return false;
07248 }
07249
07250
07251
07252 bool
07253 SparseComplexMatrix::all_elements_are_real (void) const
07254 {
07255 return mx_inline_all_real (nnz (), data ());
07256 }
07257
07258
07259
07260
07261
07262 bool
07263 SparseComplexMatrix::all_integers (double& max_val, double& min_val) const
07264 {
07265 octave_idx_type nel = nnz ();
07266
07267 if (nel == 0)
07268 return false;
07269
07270 max_val = std::real(data (0));
07271 min_val = std::real(data (0));
07272
07273 for (octave_idx_type i = 0; i < nel; i++)
07274 {
07275 Complex val = data (i);
07276
07277 double r_val = std::real (val);
07278 double i_val = std::imag (val);
07279
07280 if (r_val > max_val)
07281 max_val = r_val;
07282
07283 if (i_val > max_val)
07284 max_val = i_val;
07285
07286 if (r_val < min_val)
07287 min_val = r_val;
07288
07289 if (i_val < min_val)
07290 min_val = i_val;
07291
07292 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
07293 return false;
07294 }
07295
07296 return true;
07297 }
07298
07299 bool
07300 SparseComplexMatrix::too_large_for_float (void) const
07301 {
07302 octave_idx_type nel = nnz ();
07303
07304 for (octave_idx_type i = 0; i < nel; i++)
07305 {
07306 Complex val = data (i);
07307
07308 double r_val = std::real (val);
07309 double i_val = std::imag (val);
07310
07311 if (r_val > FLT_MAX
07312 || i_val > FLT_MAX
07313 || r_val < FLT_MIN
07314 || i_val < FLT_MIN)
07315 return true;
07316 }
07317
07318 return false;
07319 }
07320
07321
07322
07323
07324 SparseBoolMatrix
07325 SparseComplexMatrix::all (int dim) const
07326 {
07327 SPARSE_ALL_OP (dim);
07328 }
07329
07330 SparseBoolMatrix
07331 SparseComplexMatrix::any (int dim) const
07332 {
07333 SPARSE_ANY_OP (dim);
07334 }
07335
07336 SparseComplexMatrix
07337 SparseComplexMatrix::cumprod (int dim) const
07338 {
07339 SPARSE_CUMPROD (SparseComplexMatrix, Complex, cumprod);
07340 }
07341
07342 SparseComplexMatrix
07343 SparseComplexMatrix::cumsum (int dim) const
07344 {
07345 SPARSE_CUMSUM (SparseComplexMatrix, Complex, cumsum);
07346 }
07347
07348 SparseComplexMatrix
07349 SparseComplexMatrix::prod (int dim) const
07350 {
07351 if ((rows() == 1 && dim == -1) || dim == 1)
07352 return transpose (). prod (0). transpose();
07353 else
07354 {
07355 SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, *=,
07356 (cidx(j+1) - cidx(j) < nr ? 0.0 : 1.0), 1.0);
07357 }
07358 }
07359
07360 SparseComplexMatrix
07361 SparseComplexMatrix::sum (int dim) const
07362 {
07363 SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, +=, 0.0, 0.0);
07364 }
07365
07366 SparseComplexMatrix
07367 SparseComplexMatrix::sumsq (int dim) const
07368 {
07369 #define ROW_EXPR \
07370 Complex d = data (i); \
07371 tmp [ridx(i)] += d * conj (d)
07372
07373 #define COL_EXPR \
07374 Complex d = data (i); \
07375 tmp [j] += d * conj (d)
07376
07377 SPARSE_BASE_REDUCTION_OP (SparseComplexMatrix, Complex, ROW_EXPR,
07378 COL_EXPR, 0.0, 0.0);
07379
07380 #undef ROW_EXPR
07381 #undef COL_EXPR
07382 }
07383
07384 SparseMatrix SparseComplexMatrix::abs (void) const
07385 {
07386 octave_idx_type nz = nnz ();
07387 octave_idx_type nc = cols ();
07388
07389 SparseMatrix retval (rows(), nc, nz);
07390
07391 for (octave_idx_type i = 0; i < nc + 1; i++)
07392 retval.cidx (i) = cidx (i);
07393
07394 for (octave_idx_type i = 0; i < nz; i++)
07395 {
07396 retval.data (i) = std::abs (data (i));
07397 retval.ridx (i) = ridx (i);
07398 }
07399
07400 return retval;
07401 }
07402
07403 SparseComplexMatrix
07404 SparseComplexMatrix::diag (octave_idx_type k) const
07405 {
07406 return MSparse<Complex>::diag (k);
07407 }
07408
07409 std::ostream&
07410 operator << (std::ostream& os, const SparseComplexMatrix& a)
07411 {
07412 octave_idx_type nc = a.cols ();
07413
07414
07415
07416 for (octave_idx_type j = 0; j < nc; j++)
07417 {
07418 octave_quit ();
07419 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
07420 {
07421 os << a.ridx(i) + 1 << " " << j + 1 << " ";
07422 octave_write_complex (os, a.data(i));
07423 os << "\n";
07424 }
07425 }
07426
07427 return os;
07428 }
07429
07430 std::istream&
07431 operator >> (std::istream& is, SparseComplexMatrix& a)
07432 {
07433 typedef SparseComplexMatrix::element_type elt_type;
07434
07435 return read_sparse_matrix<elt_type> (is, a, octave_read_value<Complex>);
07436 }
07437
07438 SparseComplexMatrix
07439 operator * (const SparseComplexMatrix& m, const SparseMatrix& a)
07440 {
07441 SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, double);
07442 }
07443
07444 SparseComplexMatrix
07445 operator * (const SparseMatrix& m, const SparseComplexMatrix& a)
07446 {
07447 SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
07448 }
07449
07450 SparseComplexMatrix
07451 operator * (const SparseComplexMatrix& m, const SparseComplexMatrix& a)
07452 {
07453 SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
07454 }
07455
07456 ComplexMatrix
07457 operator * (const ComplexMatrix& m, const SparseMatrix& a)
07458 {
07459 FULL_SPARSE_MUL (ComplexMatrix, double, Complex (0.,0.));
07460 }
07461
07462 ComplexMatrix
07463 operator * (const Matrix& m, const SparseComplexMatrix& a)
07464 {
07465 FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
07466 }
07467
07468 ComplexMatrix
07469 operator * (const ComplexMatrix& m, const SparseComplexMatrix& a)
07470 {
07471 FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
07472 }
07473
07474 ComplexMatrix
07475 mul_trans (const ComplexMatrix& m, const SparseComplexMatrix& a)
07476 {
07477 FULL_SPARSE_MUL_TRANS (ComplexMatrix, Complex, Complex (0.,0.), );
07478 }
07479
07480 ComplexMatrix
07481 mul_herm (const ComplexMatrix& m, const SparseComplexMatrix& a)
07482 {
07483 FULL_SPARSE_MUL_TRANS (ComplexMatrix, Complex, Complex (0.,0.), conj);
07484 }
07485
07486 ComplexMatrix
07487 operator * (const SparseComplexMatrix& m, const Matrix& a)
07488 {
07489 SPARSE_FULL_MUL (ComplexMatrix, double, Complex (0.,0.));
07490 }
07491
07492 ComplexMatrix
07493 operator * (const SparseMatrix& m, const ComplexMatrix& a)
07494 {
07495 SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
07496 }
07497
07498 ComplexMatrix
07499 operator * (const SparseComplexMatrix& m, const ComplexMatrix& a)
07500 {
07501 SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
07502 }
07503
07504 ComplexMatrix
07505 trans_mul (const SparseComplexMatrix& m, const ComplexMatrix& a)
07506 {
07507 SPARSE_FULL_TRANS_MUL (ComplexMatrix, Complex, Complex (0.,0.), );
07508 }
07509
07510 ComplexMatrix
07511 herm_mul (const SparseComplexMatrix& m, const ComplexMatrix& a)
07512 {
07513 SPARSE_FULL_TRANS_MUL (ComplexMatrix, Complex, Complex (0.,0.), conj);
07514 }
07515
07516
07517 SparseComplexMatrix
07518 operator * (const DiagMatrix& d, const SparseComplexMatrix& a)
07519 {
07520 return do_mul_dm_sm<SparseComplexMatrix> (d, a);
07521 }
07522 SparseComplexMatrix
07523 operator * (const SparseComplexMatrix& a, const DiagMatrix& d)
07524 {
07525 return do_mul_sm_dm<SparseComplexMatrix> (a, d);
07526 }
07527
07528 SparseComplexMatrix
07529 operator * (const ComplexDiagMatrix& d, const SparseMatrix& a)
07530 {
07531 return do_mul_dm_sm<SparseComplexMatrix> (d, a);
07532 }
07533 SparseComplexMatrix
07534 operator * (const SparseMatrix& a, const ComplexDiagMatrix& d)
07535 {
07536 return do_mul_sm_dm<SparseComplexMatrix> (a, d);
07537 }
07538
07539 SparseComplexMatrix
07540 operator * (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
07541 {
07542 return do_mul_dm_sm<SparseComplexMatrix> (d, a);
07543 }
07544 SparseComplexMatrix
07545 operator * (const SparseComplexMatrix& a, const ComplexDiagMatrix& d)
07546 {
07547 return do_mul_sm_dm<SparseComplexMatrix> (a, d);
07548 }
07549
07550 SparseComplexMatrix
07551 operator + (const ComplexDiagMatrix& d, const SparseMatrix& a)
07552 {
07553 return do_add_dm_sm<SparseComplexMatrix> (d, a);
07554 }
07555 SparseComplexMatrix
07556 operator + (const DiagMatrix& d, const SparseComplexMatrix& a)
07557 {
07558 return do_add_dm_sm<SparseComplexMatrix> (d, a);
07559 }
07560 SparseComplexMatrix
07561 operator + (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
07562 {
07563 return do_add_dm_sm<SparseComplexMatrix> (d, a);
07564 }
07565 SparseComplexMatrix
07566 operator + (const SparseMatrix& a, const ComplexDiagMatrix& d)
07567 {
07568 return do_add_sm_dm<SparseComplexMatrix> (a, d);
07569 }
07570 SparseComplexMatrix
07571 operator + (const SparseComplexMatrix& a, const DiagMatrix& d)
07572 {
07573 return do_add_sm_dm<SparseComplexMatrix> (a, d);
07574 }
07575 SparseComplexMatrix
07576 operator + (const SparseComplexMatrix&a, const ComplexDiagMatrix& d)
07577 {
07578 return do_add_sm_dm<SparseComplexMatrix> (a, d);
07579 }
07580
07581 SparseComplexMatrix
07582 operator - (const ComplexDiagMatrix& d, const SparseMatrix& a)
07583 {
07584 return do_sub_dm_sm<SparseComplexMatrix> (d, a);
07585 }
07586 SparseComplexMatrix
07587 operator - (const DiagMatrix& d, const SparseComplexMatrix& a)
07588 {
07589 return do_sub_dm_sm<SparseComplexMatrix> (d, a);
07590 }
07591 SparseComplexMatrix
07592 operator - (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
07593 {
07594 return do_sub_dm_sm<SparseComplexMatrix> (d, a);
07595 }
07596 SparseComplexMatrix
07597 operator - (const SparseMatrix& a, const ComplexDiagMatrix& d)
07598 {
07599 return do_sub_sm_dm<SparseComplexMatrix> (a, d);
07600 }
07601 SparseComplexMatrix
07602 operator - (const SparseComplexMatrix& a, const DiagMatrix& d)
07603 {
07604 return do_sub_sm_dm<SparseComplexMatrix> (a, d);
07605 }
07606 SparseComplexMatrix
07607 operator - (const SparseComplexMatrix&a, const ComplexDiagMatrix& d)
07608 {
07609 return do_sub_sm_dm<SparseComplexMatrix> (a, d);
07610 }
07611
07612
07613
07614 SparseComplexMatrix
07615 operator * (const PermMatrix& p, const SparseComplexMatrix& a)
07616 {
07617 return octinternal_do_mul_pm_sm (p, a);
07618 }
07619
07620 SparseComplexMatrix
07621 operator * (const SparseComplexMatrix& a, const PermMatrix& p)
07622 {
07623 return octinternal_do_mul_sm_pm (a, p);
07624 }
07625
07626
07627
07628
07629 #define EMPTY_RETURN_CHECK(T) \
07630 if (nr == 0 || nc == 0) \
07631 return T (nr, nc);
07632
07633 SparseComplexMatrix
07634 min (const Complex& c, const SparseComplexMatrix& m)
07635 {
07636 SparseComplexMatrix result;
07637
07638 octave_idx_type nr = m.rows ();
07639 octave_idx_type nc = m.columns ();
07640
07641 EMPTY_RETURN_CHECK (SparseComplexMatrix);
07642
07643 if (abs(c) == 0.)
07644 return SparseComplexMatrix (nr, nc);
07645 else
07646 {
07647 result = SparseComplexMatrix (m);
07648
07649 for (octave_idx_type j = 0; j < nc; j++)
07650 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
07651 result.data(i) = xmin(c, m.data(i));
07652 }
07653
07654 return result;
07655 }
07656
07657 SparseComplexMatrix
07658 min (const SparseComplexMatrix& m, const Complex& c)
07659 {
07660 return min (c, m);
07661 }
07662
07663 SparseComplexMatrix
07664 min (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
07665 {
07666 SparseComplexMatrix r;
07667
07668 if ((a.rows() == b.rows()) && (a.cols() == b.cols()))
07669 {
07670 octave_idx_type a_nr = a.rows ();
07671 octave_idx_type a_nc = a.cols ();
07672
07673 octave_idx_type b_nr = b.rows ();
07674 octave_idx_type b_nc = b.cols ();
07675
07676 if (a_nr == 0 || b_nc == 0 || a.nnz () == 0 || b.nnz () == 0)
07677 return SparseComplexMatrix (a_nr, a_nc);
07678
07679 if (a_nr != b_nr || a_nc != b_nc)
07680 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
07681 else
07682 {
07683 r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
07684
07685 octave_idx_type jx = 0;
07686 r.cidx (0) = 0;
07687 for (octave_idx_type i = 0 ; i < a_nc ; i++)
07688 {
07689 octave_idx_type ja = a.cidx(i);
07690 octave_idx_type ja_max = a.cidx(i+1);
07691 bool ja_lt_max= ja < ja_max;
07692
07693 octave_idx_type jb = b.cidx(i);
07694 octave_idx_type jb_max = b.cidx(i+1);
07695 bool jb_lt_max = jb < jb_max;
07696
07697 while (ja_lt_max || jb_lt_max )
07698 {
07699 octave_quit ();
07700 if ((! jb_lt_max) ||
07701 (ja_lt_max && (a.ridx(ja) < b.ridx(jb))))
07702 {
07703 Complex tmp = xmin (a.data(ja), 0.);
07704 if (tmp != 0.)
07705 {
07706 r.ridx(jx) = a.ridx(ja);
07707 r.data(jx) = tmp;
07708 jx++;
07709 }
07710 ja++;
07711 ja_lt_max= ja < ja_max;
07712 }
07713 else if (( !ja_lt_max ) ||
07714 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) )
07715 {
07716 Complex tmp = xmin (0., b.data(jb));
07717 if (tmp != 0.)
07718 {
07719 r.ridx(jx) = b.ridx(jb);
07720 r.data(jx) = tmp;
07721 jx++;
07722 }
07723 jb++;
07724 jb_lt_max= jb < jb_max;
07725 }
07726 else
07727 {
07728 Complex tmp = xmin (a.data(ja), b.data(jb));
07729 if (tmp != 0.)
07730 {
07731 r.data(jx) = tmp;
07732 r.ridx(jx) = a.ridx(ja);
07733 jx++;
07734 }
07735 ja++;
07736 ja_lt_max= ja < ja_max;
07737 jb++;
07738 jb_lt_max= jb < jb_max;
07739 }
07740 }
07741 r.cidx(i+1) = jx;
07742 }
07743
07744 r.maybe_compress ();
07745 }
07746 }
07747 else
07748 (*current_liboctave_error_handler) ("matrix size mismatch");
07749
07750 return r;
07751 }
07752
07753 SparseComplexMatrix
07754 max (const Complex& c, const SparseComplexMatrix& m)
07755 {
07756 SparseComplexMatrix result;
07757
07758 octave_idx_type nr = m.rows ();
07759 octave_idx_type nc = m.columns ();
07760
07761 EMPTY_RETURN_CHECK (SparseComplexMatrix);
07762
07763
07764 if (xmax(c, 0.) != 0.)
07765 {
07766 result = SparseComplexMatrix (nr, nc, c);
07767 for (octave_idx_type j = 0; j < nc; j++)
07768 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
07769 result.xdata(m.ridx(i) + j * nr) = xmax (c, m.data(i));
07770 }
07771 else
07772 result = SparseComplexMatrix (m);
07773
07774 return result;
07775 }
07776
07777 SparseComplexMatrix
07778 max (const SparseComplexMatrix& m, const Complex& c)
07779 {
07780 return max (c, m);
07781 }
07782
07783 SparseComplexMatrix
07784 max (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
07785 {
07786 SparseComplexMatrix r;
07787
07788 if ((a.rows() == b.rows()) && (a.cols() == b.cols()))
07789 {
07790 octave_idx_type a_nr = a.rows ();
07791 octave_idx_type a_nc = a.cols ();
07792
07793 octave_idx_type b_nr = b.rows ();
07794 octave_idx_type b_nc = b.cols ();
07795
07796 if (a_nr == 0 || b_nc == 0)
07797 return SparseComplexMatrix (a_nr, a_nc);
07798 if (a.nnz () == 0)
07799 return SparseComplexMatrix (b);
07800 if (b.nnz () == 0)
07801 return SparseComplexMatrix (a);
07802
07803 if (a_nr != b_nr || a_nc != b_nc)
07804 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
07805 else
07806 {
07807 r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
07808
07809 octave_idx_type jx = 0;
07810 r.cidx (0) = 0;
07811 for (octave_idx_type i = 0 ; i < a_nc ; i++)
07812 {
07813 octave_idx_type ja = a.cidx(i);
07814 octave_idx_type ja_max = a.cidx(i+1);
07815 bool ja_lt_max= ja < ja_max;
07816
07817 octave_idx_type jb = b.cidx(i);
07818 octave_idx_type jb_max = b.cidx(i+1);
07819 bool jb_lt_max = jb < jb_max;
07820
07821 while (ja_lt_max || jb_lt_max )
07822 {
07823 octave_quit ();
07824 if ((! jb_lt_max) ||
07825 (ja_lt_max && (a.ridx(ja) < b.ridx(jb))))
07826 {
07827 Complex tmp = xmax (a.data(ja), 0.);
07828 if (tmp != 0.)
07829 {
07830 r.ridx(jx) = a.ridx(ja);
07831 r.data(jx) = tmp;
07832 jx++;
07833 }
07834 ja++;
07835 ja_lt_max= ja < ja_max;
07836 }
07837 else if (( !ja_lt_max ) ||
07838 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) )
07839 {
07840 Complex tmp = xmax (0., b.data(jb));
07841 if (tmp != 0.)
07842 {
07843 r.ridx(jx) = b.ridx(jb);
07844 r.data(jx) = tmp;
07845 jx++;
07846 }
07847 jb++;
07848 jb_lt_max= jb < jb_max;
07849 }
07850 else
07851 {
07852 Complex tmp = xmax (a.data(ja), b.data(jb));
07853 if (tmp != 0.)
07854 {
07855 r.data(jx) = tmp;
07856 r.ridx(jx) = a.ridx(ja);
07857 jx++;
07858 }
07859 ja++;
07860 ja_lt_max= ja < ja_max;
07861 jb++;
07862 jb_lt_max= jb < jb_max;
07863 }
07864 }
07865 r.cidx(i+1) = jx;
07866 }
07867
07868 r.maybe_compress ();
07869 }
07870 }
07871 else
07872 (*current_liboctave_error_handler) ("matrix size mismatch");
07873
07874 return r;
07875 }
07876
07877 SPARSE_SMS_CMP_OPS (SparseComplexMatrix, 0.0, real, Complex,
07878 0.0, real)
07879 SPARSE_SMS_BOOL_OPS (SparseComplexMatrix, Complex, 0.0)
07880
07881 SPARSE_SSM_CMP_OPS (Complex, 0.0, real, SparseComplexMatrix,
07882 0.0, real)
07883 SPARSE_SSM_BOOL_OPS (Complex, SparseComplexMatrix, 0.0)
07884
07885 SPARSE_SMSM_CMP_OPS (SparseComplexMatrix, 0.0, real, SparseComplexMatrix,
07886 0.0, real)
07887 SPARSE_SMSM_BOOL_OPS (SparseComplexMatrix, SparseComplexMatrix, 0.0)
1.7.1
