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