[9434] | 1 | /// \ingroup newmat
|
---|
| 2 | ///@{
|
---|
| 3 |
|
---|
| 4 | /// \file newmat6.cpp
|
---|
| 5 | /// Operators, element access.
|
---|
| 6 |
|
---|
| 7 | // Copyright (C) 1991,2,3,4: R B Davies
|
---|
| 8 |
|
---|
| 9 | #include "include.h"
|
---|
| 10 |
|
---|
| 11 | #include "newmat.h"
|
---|
| 12 | #include "newmatrc.h"
|
---|
| 13 |
|
---|
| 14 | #ifdef use_namespace
|
---|
| 15 | namespace NEWMAT {
|
---|
| 16 | #endif
|
---|
| 17 |
|
---|
| 18 |
|
---|
| 19 |
|
---|
| 20 | #ifdef DO_REPORT
|
---|
| 21 | #define REPORT { static ExeCounter ExeCount(__LINE__,6); ++ExeCount; }
|
---|
| 22 | #else
|
---|
| 23 | #define REPORT {}
|
---|
| 24 | #endif
|
---|
| 25 |
|
---|
| 26 | /*************************** general utilities *************************/
|
---|
| 27 |
|
---|
| 28 | static int tristore(int n) // els in triangular matrix
|
---|
| 29 | { return (n*(n+1))/2; }
|
---|
| 30 |
|
---|
| 31 |
|
---|
| 32 | /****************************** operators *******************************/
|
---|
| 33 |
|
---|
| 34 | Real& Matrix::operator()(int m, int n)
|
---|
| 35 | {
|
---|
| 36 | REPORT
|
---|
| 37 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val)
|
---|
| 38 | Throw(IndexException(m,n,*this));
|
---|
| 39 | return store[(m-1)*ncols_val+n-1];
|
---|
| 40 | }
|
---|
| 41 |
|
---|
| 42 | Real& SymmetricMatrix::operator()(int m, int n)
|
---|
| 43 | {
|
---|
| 44 | REPORT
|
---|
| 45 | if (m<=0 || n<=0 || m>nrows_val || n>ncols_val)
|
---|
| 46 | Throw(IndexException(m,n,*this));
|
---|
| 47 | if (m>=n) return store[tristore(m-1)+n-1];
|
---|
| 48 | else return store[tristore(n-1)+m-1];
|
---|
| 49 | }
|
---|
| 50 |
|
---|
| 51 | Real& UpperTriangularMatrix::operator()(int m, int n)
|
---|
| 52 | {
|
---|
| 53 | REPORT
|
---|
| 54 | if (m<=0 || n<m || n>ncols_val)
|
---|
| 55 | Throw(IndexException(m,n,*this));
|
---|
| 56 | return store[(m-1)*ncols_val+n-1-tristore(m-1)];
|
---|
| 57 | }
|
---|
| 58 |
|
---|
| 59 | Real& LowerTriangularMatrix::operator()(int m, int n)
|
---|
| 60 | {
|
---|
| 61 | REPORT
|
---|
| 62 | if (n<=0 || m<n || m>nrows_val)
|
---|
| 63 | Throw(IndexException(m,n,*this));
|
---|
| 64 | return store[tristore(m-1)+n-1];
|
---|
| 65 | }
|
---|
| 66 |
|
---|
| 67 | Real& DiagonalMatrix::operator()(int m, int n)
|
---|
| 68 | {
|
---|
| 69 | REPORT
|
---|
| 70 | if (n<=0 || m!=n || m>nrows_val || n>ncols_val)
|
---|
| 71 | Throw(IndexException(m,n,*this));
|
---|
| 72 | return store[n-1];
|
---|
| 73 | }
|
---|
| 74 |
|
---|
| 75 | Real& DiagonalMatrix::operator()(int m)
|
---|
| 76 | {
|
---|
| 77 | REPORT
|
---|
| 78 | if (m<=0 || m>nrows_val) Throw(IndexException(m,*this));
|
---|
| 79 | return store[m-1];
|
---|
| 80 | }
|
---|
| 81 |
|
---|
| 82 | Real& ColumnVector::operator()(int m)
|
---|
| 83 | {
|
---|
| 84 | REPORT
|
---|
| 85 | if (m<=0 || m> nrows_val) Throw(IndexException(m,*this));
|
---|
| 86 | return store[m-1];
|
---|
| 87 | }
|
---|
| 88 |
|
---|
| 89 | Real& RowVector::operator()(int n)
|
---|
| 90 | {
|
---|
| 91 | REPORT
|
---|
| 92 | if (n<=0 || n> ncols_val) Throw(IndexException(n,*this));
|
---|
| 93 | return store[n-1];
|
---|
| 94 | }
|
---|
| 95 |
|
---|
| 96 | Real& BandMatrix::operator()(int m, int n)
|
---|
| 97 | {
|
---|
| 98 | REPORT
|
---|
| 99 | int w = upper_val+lower_val+1; int i = lower_val+n-m;
|
---|
| 100 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
|
---|
| 101 | Throw(IndexException(m,n,*this));
|
---|
| 102 | return store[w*(m-1)+i];
|
---|
| 103 | }
|
---|
| 104 |
|
---|
| 105 | Real& UpperBandMatrix::operator()(int m, int n)
|
---|
| 106 | {
|
---|
| 107 | REPORT
|
---|
| 108 | int w = upper_val+1; int i = n-m;
|
---|
| 109 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
|
---|
| 110 | Throw(IndexException(m,n,*this));
|
---|
| 111 | return store[w*(m-1)+i];
|
---|
| 112 | }
|
---|
| 113 |
|
---|
| 114 | Real& LowerBandMatrix::operator()(int m, int n)
|
---|
| 115 | {
|
---|
| 116 | REPORT
|
---|
| 117 | int w = lower_val+1; int i = lower_val+n-m;
|
---|
| 118 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
|
---|
| 119 | Throw(IndexException(m,n,*this));
|
---|
| 120 | return store[w*(m-1)+i];
|
---|
| 121 | }
|
---|
| 122 |
|
---|
| 123 | Real& SymmetricBandMatrix::operator()(int m, int n)
|
---|
| 124 | {
|
---|
| 125 | REPORT
|
---|
| 126 | int w = lower_val+1;
|
---|
| 127 | if (m>=n)
|
---|
| 128 | {
|
---|
| 129 | REPORT
|
---|
| 130 | int i = lower_val+n-m;
|
---|
| 131 | if ( m>nrows_val || n<=0 || i<0 )
|
---|
| 132 | Throw(IndexException(m,n,*this));
|
---|
| 133 | return store[w*(m-1)+i];
|
---|
| 134 | }
|
---|
| 135 | else
|
---|
| 136 | {
|
---|
| 137 | REPORT
|
---|
| 138 | int i = lower_val+m-n;
|
---|
| 139 | if ( n>nrows_val || m<=0 || i<0 )
|
---|
| 140 | Throw(IndexException(m,n,*this));
|
---|
| 141 | return store[w*(n-1)+i];
|
---|
| 142 | }
|
---|
| 143 | }
|
---|
| 144 |
|
---|
| 145 |
|
---|
| 146 | Real Matrix::operator()(int m, int n) const
|
---|
| 147 | {
|
---|
| 148 | REPORT
|
---|
| 149 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val)
|
---|
| 150 | Throw(IndexException(m,n,*this));
|
---|
| 151 | return store[(m-1)*ncols_val+n-1];
|
---|
| 152 | }
|
---|
| 153 |
|
---|
| 154 | Real SymmetricMatrix::operator()(int m, int n) const
|
---|
| 155 | {
|
---|
| 156 | REPORT
|
---|
| 157 | if (m<=0 || n<=0 || m>nrows_val || n>ncols_val)
|
---|
| 158 | Throw(IndexException(m,n,*this));
|
---|
| 159 | if (m>=n) return store[tristore(m-1)+n-1];
|
---|
| 160 | else return store[tristore(n-1)+m-1];
|
---|
| 161 | }
|
---|
| 162 |
|
---|
| 163 | Real UpperTriangularMatrix::operator()(int m, int n) const
|
---|
| 164 | {
|
---|
| 165 | REPORT
|
---|
| 166 | if (m<=0 || n<m || n>ncols_val)
|
---|
| 167 | Throw(IndexException(m,n,*this));
|
---|
| 168 | return store[(m-1)*ncols_val+n-1-tristore(m-1)];
|
---|
| 169 | }
|
---|
| 170 |
|
---|
| 171 | Real LowerTriangularMatrix::operator()(int m, int n) const
|
---|
| 172 | {
|
---|
| 173 | REPORT
|
---|
| 174 | if (n<=0 || m<n || m>nrows_val)
|
---|
| 175 | Throw(IndexException(m,n,*this));
|
---|
| 176 | return store[tristore(m-1)+n-1];
|
---|
| 177 | }
|
---|
| 178 |
|
---|
| 179 | Real DiagonalMatrix::operator()(int m, int n) const
|
---|
| 180 | {
|
---|
| 181 | REPORT
|
---|
| 182 | if (n<=0 || m!=n || m>nrows_val || n>ncols_val)
|
---|
| 183 | Throw(IndexException(m,n,*this));
|
---|
| 184 | return store[n-1];
|
---|
| 185 | }
|
---|
| 186 |
|
---|
| 187 | Real DiagonalMatrix::operator()(int m) const
|
---|
| 188 | {
|
---|
| 189 | REPORT
|
---|
| 190 | if (m<=0 || m>nrows_val) Throw(IndexException(m,*this));
|
---|
| 191 | return store[m-1];
|
---|
| 192 | }
|
---|
| 193 |
|
---|
| 194 | Real ColumnVector::operator()(int m) const
|
---|
| 195 | {
|
---|
| 196 | REPORT
|
---|
| 197 | if (m<=0 || m> nrows_val) Throw(IndexException(m,*this));
|
---|
| 198 | return store[m-1];
|
---|
| 199 | }
|
---|
| 200 |
|
---|
| 201 | Real RowVector::operator()(int n) const
|
---|
| 202 | {
|
---|
| 203 | REPORT
|
---|
| 204 | if (n<=0 || n> ncols_val) Throw(IndexException(n,*this));
|
---|
| 205 | return store[n-1];
|
---|
| 206 | }
|
---|
| 207 |
|
---|
| 208 | Real BandMatrix::operator()(int m, int n) const
|
---|
| 209 | {
|
---|
| 210 | REPORT
|
---|
| 211 | int w = upper_val+lower_val+1; int i = lower_val+n-m;
|
---|
| 212 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
|
---|
| 213 | Throw(IndexException(m,n,*this));
|
---|
| 214 | return store[w*(m-1)+i];
|
---|
| 215 | }
|
---|
| 216 |
|
---|
| 217 | Real UpperBandMatrix::operator()(int m, int n) const
|
---|
| 218 | {
|
---|
| 219 | REPORT
|
---|
| 220 | int w = upper_val+1; int i = n-m;
|
---|
| 221 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
|
---|
| 222 | Throw(IndexException(m,n,*this));
|
---|
| 223 | return store[w*(m-1)+i];
|
---|
| 224 | }
|
---|
| 225 |
|
---|
| 226 | Real LowerBandMatrix::operator()(int m, int n) const
|
---|
| 227 | {
|
---|
| 228 | REPORT
|
---|
| 229 | int w = lower_val+1; int i = lower_val+n-m;
|
---|
| 230 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
|
---|
| 231 | Throw(IndexException(m,n,*this));
|
---|
| 232 | return store[w*(m-1)+i];
|
---|
| 233 | }
|
---|
| 234 |
|
---|
| 235 | Real SymmetricBandMatrix::operator()(int m, int n) const
|
---|
| 236 | {
|
---|
| 237 | REPORT
|
---|
| 238 | int w = lower_val+1;
|
---|
| 239 | if (m>=n)
|
---|
| 240 | {
|
---|
| 241 | REPORT
|
---|
| 242 | int i = lower_val+n-m;
|
---|
| 243 | if ( m>nrows_val || n<=0 || i<0 )
|
---|
| 244 | Throw(IndexException(m,n,*this));
|
---|
| 245 | return store[w*(m-1)+i];
|
---|
| 246 | }
|
---|
| 247 | else
|
---|
| 248 | {
|
---|
| 249 | REPORT
|
---|
| 250 | int i = lower_val+m-n;
|
---|
| 251 | if ( n>nrows_val || m<=0 || i<0 )
|
---|
| 252 | Throw(IndexException(m,n,*this));
|
---|
| 253 | return store[w*(n-1)+i];
|
---|
| 254 | }
|
---|
| 255 | }
|
---|
| 256 |
|
---|
| 257 |
|
---|
| 258 | Real BaseMatrix::as_scalar() const
|
---|
| 259 | {
|
---|
| 260 | REPORT
|
---|
| 261 | GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
|
---|
| 262 |
|
---|
| 263 | if (gm->nrows_val!=1 || gm->ncols_val!=1)
|
---|
| 264 | {
|
---|
| 265 | Tracer tr("as_scalar");
|
---|
| 266 | Try
|
---|
| 267 | { Throw(ProgramException("Cannot convert to scalar", *gm)); }
|
---|
| 268 | CatchAll { gm->tDelete(); ReThrow; }
|
---|
| 269 | }
|
---|
| 270 |
|
---|
| 271 | Real x = *(gm->store); gm->tDelete(); return x;
|
---|
| 272 | }
|
---|
| 273 |
|
---|
| 274 |
|
---|
| 275 | AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
|
---|
| 276 | { REPORT return AddedMatrix(this, &bm); }
|
---|
| 277 |
|
---|
| 278 | SPMatrix SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
|
---|
| 279 | { REPORT return SPMatrix(&bm1, &bm2); }
|
---|
| 280 |
|
---|
| 281 | KPMatrix KP(const BaseMatrix& bm1,const BaseMatrix& bm2)
|
---|
| 282 | { REPORT return KPMatrix(&bm1, &bm2); }
|
---|
| 283 |
|
---|
| 284 | MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
|
---|
| 285 | { REPORT return MultipliedMatrix(this, &bm); }
|
---|
| 286 |
|
---|
| 287 | ConcatenatedMatrix BaseMatrix::operator|(const BaseMatrix& bm) const
|
---|
| 288 | { REPORT return ConcatenatedMatrix(this, &bm); }
|
---|
| 289 |
|
---|
| 290 | StackedMatrix BaseMatrix::operator&(const BaseMatrix& bm) const
|
---|
| 291 | { REPORT return StackedMatrix(this, &bm); }
|
---|
| 292 |
|
---|
| 293 | SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
|
---|
| 294 | { REPORT return SolvedMatrix(bm, &bmx); }
|
---|
| 295 |
|
---|
| 296 | SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
|
---|
| 297 | { REPORT return SubtractedMatrix(this, &bm); }
|
---|
| 298 |
|
---|
| 299 | ShiftedMatrix BaseMatrix::operator+(Real f) const
|
---|
| 300 | { REPORT return ShiftedMatrix(this, f); }
|
---|
| 301 |
|
---|
| 302 | ShiftedMatrix operator+(Real f, const BaseMatrix& BM)
|
---|
| 303 | { REPORT return ShiftedMatrix(&BM, f); }
|
---|
| 304 |
|
---|
| 305 | NegShiftedMatrix operator-(Real f, const BaseMatrix& bm)
|
---|
| 306 | { REPORT return NegShiftedMatrix(f, &bm); }
|
---|
| 307 |
|
---|
| 308 | ScaledMatrix BaseMatrix::operator*(Real f) const
|
---|
| 309 | { REPORT return ScaledMatrix(this, f); }
|
---|
| 310 |
|
---|
| 311 | ScaledMatrix BaseMatrix::operator/(Real f) const
|
---|
| 312 | { REPORT return ScaledMatrix(this, 1.0/f); }
|
---|
| 313 |
|
---|
| 314 | ScaledMatrix operator*(Real f, const BaseMatrix& BM)
|
---|
| 315 | { REPORT return ScaledMatrix(&BM, f); }
|
---|
| 316 |
|
---|
| 317 | ShiftedMatrix BaseMatrix::operator-(Real f) const
|
---|
| 318 | { REPORT return ShiftedMatrix(this, -f); }
|
---|
| 319 |
|
---|
| 320 | TransposedMatrix BaseMatrix::t() const
|
---|
| 321 | { REPORT return TransposedMatrix(this); }
|
---|
| 322 |
|
---|
| 323 | NegatedMatrix BaseMatrix::operator-() const
|
---|
| 324 | { REPORT return NegatedMatrix(this); }
|
---|
| 325 |
|
---|
| 326 | ReversedMatrix BaseMatrix::reverse() const
|
---|
| 327 | { REPORT return ReversedMatrix(this); }
|
---|
| 328 |
|
---|
| 329 | InvertedMatrix BaseMatrix::i() const
|
---|
| 330 | { REPORT return InvertedMatrix(this); }
|
---|
| 331 |
|
---|
| 332 |
|
---|
| 333 | RowedMatrix BaseMatrix::as_row() const
|
---|
| 334 | { REPORT return RowedMatrix(this); }
|
---|
| 335 |
|
---|
| 336 | ColedMatrix BaseMatrix::as_column() const
|
---|
| 337 | { REPORT return ColedMatrix(this); }
|
---|
| 338 |
|
---|
| 339 | DiagedMatrix BaseMatrix::as_diagonal() const
|
---|
| 340 | { REPORT return DiagedMatrix(this); }
|
---|
| 341 |
|
---|
| 342 | MatedMatrix BaseMatrix::as_matrix(int nrx, int ncx) const
|
---|
| 343 | { REPORT return MatedMatrix(this,nrx,ncx); }
|
---|
| 344 |
|
---|
| 345 |
|
---|
| 346 | void GeneralMatrix::operator=(Real f)
|
---|
| 347 | { REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }
|
---|
| 348 |
|
---|
| 349 | void Matrix::operator=(const BaseMatrix& X)
|
---|
| 350 | {
|
---|
| 351 | REPORT //CheckConversion(X);
|
---|
| 352 | // MatrixConversionCheck mcc;
|
---|
| 353 | Eq(X,MatrixType::Rt);
|
---|
[9470] | 354 | }
|
---|
[9434] | 355 |
|
---|
| 356 | void SquareMatrix::operator=(const BaseMatrix& X)
|
---|
| 357 | {
|
---|
| 358 | REPORT //CheckConversion(X);
|
---|
| 359 | // MatrixConversionCheck mcc;
|
---|
| 360 | Eq(X,MatrixType::Rt);
|
---|
| 361 | if (nrows_val != ncols_val)
|
---|
| 362 | { Tracer tr("SquareMatrix(=)"); Throw(NotSquareException(*this)); }
|
---|
| 363 | }
|
---|
| 364 |
|
---|
| 365 | void SquareMatrix::operator=(const Matrix& m)
|
---|
| 366 | {
|
---|
| 367 | REPORT
|
---|
| 368 | if (m.nrows_val != m.ncols_val)
|
---|
| 369 | { Tracer tr("SquareMatrix(=Matrix)"); Throw(NotSquareException(*this)); }
|
---|
| 370 | Eq(m);
|
---|
| 371 | }
|
---|
| 372 |
|
---|
| 373 | void RowVector::operator=(const BaseMatrix& X)
|
---|
| 374 | {
|
---|
| 375 | REPORT // CheckConversion(X);
|
---|
| 376 | // MatrixConversionCheck mcc;
|
---|
| 377 | Eq(X,MatrixType::RV);
|
---|
| 378 | if (nrows_val!=1)
|
---|
| 379 | { Tracer tr("RowVector(=)"); Throw(VectorException(*this)); }
|
---|
| 380 | }
|
---|
| 381 |
|
---|
| 382 | void ColumnVector::operator=(const BaseMatrix& X)
|
---|
| 383 | {
|
---|
| 384 | REPORT //CheckConversion(X);
|
---|
| 385 | // MatrixConversionCheck mcc;
|
---|
| 386 | Eq(X,MatrixType::CV);
|
---|
| 387 | if (ncols_val!=1)
|
---|
| 388 | { Tracer tr("ColumnVector(=)"); Throw(VectorException(*this)); }
|
---|
| 389 | }
|
---|
| 390 |
|
---|
| 391 | void SymmetricMatrix::operator=(const BaseMatrix& X)
|
---|
| 392 | {
|
---|
| 393 | REPORT // CheckConversion(X);
|
---|
| 394 | // MatrixConversionCheck mcc;
|
---|
| 395 | Eq(X,MatrixType::Sm);
|
---|
| 396 | }
|
---|
| 397 |
|
---|
| 398 | void UpperTriangularMatrix::operator=(const BaseMatrix& X)
|
---|
| 399 | {
|
---|
| 400 | REPORT //CheckConversion(X);
|
---|
| 401 | // MatrixConversionCheck mcc;
|
---|
| 402 | Eq(X,MatrixType::UT);
|
---|
| 403 | }
|
---|
| 404 |
|
---|
| 405 | void LowerTriangularMatrix::operator=(const BaseMatrix& X)
|
---|
| 406 | {
|
---|
| 407 | REPORT //CheckConversion(X);
|
---|
| 408 | // MatrixConversionCheck mcc;
|
---|
| 409 | Eq(X,MatrixType::LT);
|
---|
| 410 | }
|
---|
| 411 |
|
---|
| 412 | void DiagonalMatrix::operator=(const BaseMatrix& X)
|
---|
| 413 | {
|
---|
| 414 | REPORT // CheckConversion(X);
|
---|
| 415 | // MatrixConversionCheck mcc;
|
---|
| 416 | Eq(X,MatrixType::Dg);
|
---|
| 417 | }
|
---|
| 418 |
|
---|
| 419 | void IdentityMatrix::operator=(const BaseMatrix& X)
|
---|
| 420 | {
|
---|
| 421 | REPORT // CheckConversion(X);
|
---|
| 422 | // MatrixConversionCheck mcc;
|
---|
| 423 | Eq(X,MatrixType::Id);
|
---|
| 424 | }
|
---|
| 425 |
|
---|
| 426 |
|
---|
| 427 | void CroutMatrix::operator=(const CroutMatrix& gm)
|
---|
| 428 | {
|
---|
| 429 | if (&gm == this) { REPORT tag_val = -1; return; }
|
---|
| 430 | REPORT
|
---|
[9470] | 431 | if (indx) { delete [] indx; indx = 0; }
|
---|
[9434] | 432 | ((CroutMatrix&)gm).get_aux(*this);
|
---|
| 433 | Eq(gm);
|
---|
| 434 | }
|
---|
| 435 |
|
---|
| 436 |
|
---|
| 437 |
|
---|
| 438 |
|
---|
[9470] | 439 |
|
---|
[9434] | 440 | void GeneralMatrix::operator<<(const double* r)
|
---|
| 441 | {
|
---|
| 442 | REPORT
|
---|
| 443 | int i = storage; Real* s=store;
|
---|
| 444 | while(i--) *s++ = (Real)*r++;
|
---|
| 445 | }
|
---|
| 446 |
|
---|
| 447 |
|
---|
| 448 | void GeneralMatrix::operator<<(const float* r)
|
---|
| 449 | {
|
---|
| 450 | REPORT
|
---|
| 451 | int i = storage; Real* s=store;
|
---|
| 452 | while(i--) *s++ = (Real)*r++;
|
---|
| 453 | }
|
---|
| 454 |
|
---|
| 455 |
|
---|
| 456 | void GeneralMatrix::operator<<(const int* r)
|
---|
| 457 | {
|
---|
| 458 | REPORT
|
---|
| 459 | int i = storage; Real* s=store;
|
---|
| 460 | while(i--) *s++ = (Real)*r++;
|
---|
| 461 | }
|
---|
| 462 |
|
---|
| 463 |
|
---|
| 464 | void GenericMatrix::operator=(const GenericMatrix& bmx)
|
---|
| 465 | {
|
---|
| 466 | if (&bmx != this) { REPORT if (gm) delete gm; gm = bmx.gm->Image();}
|
---|
| 467 | else { REPORT }
|
---|
| 468 | gm->Protect();
|
---|
| 469 | }
|
---|
| 470 |
|
---|
| 471 | void GenericMatrix::operator=(const BaseMatrix& bmx)
|
---|
| 472 | {
|
---|
| 473 | if (gm)
|
---|
| 474 | {
|
---|
| 475 | int counter=bmx.search(gm);
|
---|
| 476 | if (counter==0) { REPORT delete gm; gm=0; }
|
---|
| 477 | else { REPORT gm->Release(counter); }
|
---|
| 478 | }
|
---|
| 479 | else { REPORT }
|
---|
| 480 | GeneralMatrix* gmx = ((BaseMatrix&)bmx).Evaluate();
|
---|
| 481 | if (gmx != gm) { REPORT if (gm) delete gm; gm = gmx->Image(); }
|
---|
| 482 | else { REPORT }
|
---|
| 483 | gm->Protect();
|
---|
| 484 | }
|
---|
| 485 |
|
---|
| 486 |
|
---|
| 487 | /*************************** += etc ***************************************/
|
---|
| 488 |
|
---|
| 489 |
|
---|
| 490 | // GeneralMatrix operators
|
---|
| 491 |
|
---|
| 492 | void GeneralMatrix::operator+=(const BaseMatrix& X)
|
---|
| 493 | {
|
---|
| 494 | REPORT
|
---|
| 495 | Tracer tr("GeneralMatrix::operator+=");
|
---|
| 496 | // MatrixConversionCheck mcc;
|
---|
| 497 | Protect(); // so it cannot get deleted
|
---|
| 498 | // during Evaluate
|
---|
| 499 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
|
---|
| 500 | AddedMatrix am(this,gm);
|
---|
| 501 | if (gm==this) Release(2); else Release();
|
---|
| 502 | Eq2(am,type());
|
---|
| 503 | }
|
---|
| 504 |
|
---|
| 505 | void GeneralMatrix::operator-=(const BaseMatrix& X)
|
---|
| 506 | {
|
---|
| 507 | REPORT
|
---|
| 508 | Tracer tr("GeneralMatrix::operator-=");
|
---|
| 509 | // MatrixConversionCheck mcc;
|
---|
| 510 | Protect(); // so it cannot get deleted
|
---|
| 511 | // during Evaluate
|
---|
| 512 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
|
---|
| 513 | SubtractedMatrix am(this,gm);
|
---|
| 514 | if (gm==this) Release(2); else Release();
|
---|
| 515 | Eq2(am,type());
|
---|
| 516 | }
|
---|
| 517 |
|
---|
| 518 | void GeneralMatrix::operator*=(const BaseMatrix& X)
|
---|
| 519 | {
|
---|
| 520 | REPORT
|
---|
| 521 | Tracer tr("GeneralMatrix::operator*=");
|
---|
| 522 | // MatrixConversionCheck mcc;
|
---|
| 523 | Protect(); // so it cannot get deleted
|
---|
| 524 | // during Evaluate
|
---|
| 525 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
|
---|
| 526 | MultipliedMatrix am(this,gm);
|
---|
| 527 | if (gm==this) Release(2); else Release();
|
---|
| 528 | Eq2(am,type());
|
---|
| 529 | }
|
---|
| 530 |
|
---|
| 531 | void GeneralMatrix::operator|=(const BaseMatrix& X)
|
---|
| 532 | {
|
---|
| 533 | REPORT
|
---|
| 534 | Tracer tr("GeneralMatrix::operator|=");
|
---|
| 535 | // MatrixConversionCheck mcc;
|
---|
| 536 | Protect(); // so it cannot get deleted
|
---|
| 537 | // during Evaluate
|
---|
| 538 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
|
---|
| 539 | ConcatenatedMatrix am(this,gm);
|
---|
| 540 | if (gm==this) Release(2); else Release();
|
---|
| 541 | Eq2(am,type());
|
---|
| 542 | }
|
---|
| 543 |
|
---|
| 544 | void GeneralMatrix::operator&=(const BaseMatrix& X)
|
---|
| 545 | {
|
---|
| 546 | REPORT
|
---|
| 547 | Tracer tr("GeneralMatrix::operator&=");
|
---|
| 548 | // MatrixConversionCheck mcc;
|
---|
| 549 | Protect(); // so it cannot get deleted
|
---|
| 550 | // during Evaluate
|
---|
| 551 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
|
---|
| 552 | StackedMatrix am(this,gm);
|
---|
| 553 | if (gm==this) Release(2); else Release();
|
---|
| 554 | Eq2(am,type());
|
---|
| 555 | }
|
---|
| 556 |
|
---|
| 557 | void GeneralMatrix::operator+=(Real r)
|
---|
| 558 | {
|
---|
| 559 | REPORT
|
---|
| 560 | Tracer tr("GeneralMatrix::operator+=(Real)");
|
---|
| 561 | // MatrixConversionCheck mcc;
|
---|
| 562 | ShiftedMatrix am(this,r);
|
---|
| 563 | Release(); Eq2(am,type());
|
---|
| 564 | }
|
---|
| 565 |
|
---|
| 566 | void GeneralMatrix::operator*=(Real r)
|
---|
| 567 | {
|
---|
| 568 | REPORT
|
---|
| 569 | Tracer tr("GeneralMatrix::operator*=(Real)");
|
---|
| 570 | // MatrixConversionCheck mcc;
|
---|
| 571 | ScaledMatrix am(this,r);
|
---|
| 572 | Release(); Eq2(am,type());
|
---|
| 573 | }
|
---|
| 574 |
|
---|
| 575 |
|
---|
| 576 | // Generic matrix operators
|
---|
| 577 |
|
---|
| 578 | void GenericMatrix::operator+=(const BaseMatrix& X)
|
---|
| 579 | {
|
---|
| 580 | REPORT
|
---|
| 581 | Tracer tr("GenericMatrix::operator+=");
|
---|
| 582 | if (!gm) Throw(ProgramException("GenericMatrix is null"));
|
---|
| 583 | gm->Protect(); // so it cannot get deleted during Evaluate
|
---|
| 584 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
|
---|
| 585 | AddedMatrix am(gm,gmx);
|
---|
| 586 | if (gmx==gm) gm->Release(2); else gm->Release();
|
---|
| 587 | GeneralMatrix* gmy = am.Evaluate();
|
---|
| 588 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
|
---|
| 589 | else { REPORT }
|
---|
| 590 | gm->Protect();
|
---|
| 591 | }
|
---|
| 592 |
|
---|
| 593 | void GenericMatrix::operator-=(const BaseMatrix& X)
|
---|
| 594 | {
|
---|
| 595 | REPORT
|
---|
| 596 | Tracer tr("GenericMatrix::operator-=");
|
---|
| 597 | if (!gm) Throw(ProgramException("GenericMatrix is null"));
|
---|
| 598 | gm->Protect(); // so it cannot get deleted during Evaluate
|
---|
| 599 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
|
---|
| 600 | SubtractedMatrix am(gm,gmx);
|
---|
| 601 | if (gmx==gm) gm->Release(2); else gm->Release();
|
---|
| 602 | GeneralMatrix* gmy = am.Evaluate();
|
---|
| 603 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
|
---|
| 604 | else { REPORT }
|
---|
| 605 | gm->Protect();
|
---|
| 606 | }
|
---|
| 607 |
|
---|
| 608 | void GenericMatrix::operator*=(const BaseMatrix& X)
|
---|
| 609 | {
|
---|
| 610 | REPORT
|
---|
| 611 | Tracer tr("GenericMatrix::operator*=");
|
---|
| 612 | if (!gm) Throw(ProgramException("GenericMatrix is null"));
|
---|
| 613 | gm->Protect(); // so it cannot get deleted during Evaluate
|
---|
| 614 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
|
---|
| 615 | MultipliedMatrix am(gm,gmx);
|
---|
| 616 | if (gmx==gm) gm->Release(2); else gm->Release();
|
---|
| 617 | GeneralMatrix* gmy = am.Evaluate();
|
---|
| 618 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
|
---|
| 619 | else { REPORT }
|
---|
| 620 | gm->Protect();
|
---|
| 621 | }
|
---|
| 622 |
|
---|
| 623 | void GenericMatrix::operator|=(const BaseMatrix& X)
|
---|
| 624 | {
|
---|
| 625 | REPORT
|
---|
| 626 | Tracer tr("GenericMatrix::operator|=");
|
---|
| 627 | if (!gm) Throw(ProgramException("GenericMatrix is null"));
|
---|
| 628 | gm->Protect(); // so it cannot get deleted during Evaluate
|
---|
| 629 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
|
---|
| 630 | ConcatenatedMatrix am(gm,gmx);
|
---|
| 631 | if (gmx==gm) gm->Release(2); else gm->Release();
|
---|
| 632 | GeneralMatrix* gmy = am.Evaluate();
|
---|
| 633 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
|
---|
| 634 | else { REPORT }
|
---|
| 635 | gm->Protect();
|
---|
| 636 | }
|
---|
| 637 |
|
---|
| 638 | void GenericMatrix::operator&=(const BaseMatrix& X)
|
---|
| 639 | {
|
---|
| 640 | REPORT
|
---|
| 641 | Tracer tr("GenericMatrix::operator&=");
|
---|
| 642 | if (!gm) Throw(ProgramException("GenericMatrix is null"));
|
---|
| 643 | gm->Protect(); // so it cannot get deleted during Evaluate
|
---|
| 644 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
|
---|
| 645 | StackedMatrix am(gm,gmx);
|
---|
| 646 | if (gmx==gm) gm->Release(2); else gm->Release();
|
---|
| 647 | GeneralMatrix* gmy = am.Evaluate();
|
---|
| 648 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
|
---|
| 649 | else { REPORT }
|
---|
| 650 | gm->Protect();
|
---|
| 651 | }
|
---|
| 652 |
|
---|
| 653 | void GenericMatrix::operator+=(Real r)
|
---|
| 654 | {
|
---|
| 655 | REPORT
|
---|
| 656 | Tracer tr("GenericMatrix::operator+= (Real)");
|
---|
| 657 | if (!gm) Throw(ProgramException("GenericMatrix is null"));
|
---|
| 658 | ShiftedMatrix am(gm,r);
|
---|
| 659 | gm->Release();
|
---|
| 660 | GeneralMatrix* gmy = am.Evaluate();
|
---|
| 661 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
|
---|
| 662 | else { REPORT }
|
---|
| 663 | gm->Protect();
|
---|
| 664 | }
|
---|
| 665 |
|
---|
| 666 | void GenericMatrix::operator*=(Real r)
|
---|
| 667 | {
|
---|
| 668 | REPORT
|
---|
| 669 | Tracer tr("GenericMatrix::operator*= (Real)");
|
---|
| 670 | if (!gm) Throw(ProgramException("GenericMatrix is null"));
|
---|
| 671 | ScaledMatrix am(gm,r);
|
---|
| 672 | gm->Release();
|
---|
| 673 | GeneralMatrix* gmy = am.Evaluate();
|
---|
| 674 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
|
---|
| 675 | else { REPORT }
|
---|
| 676 | gm->Protect();
|
---|
| 677 | }
|
---|
| 678 |
|
---|
| 679 |
|
---|
| 680 | /************************* element access *********************************/
|
---|
| 681 |
|
---|
| 682 | Real& Matrix::element(int m, int n)
|
---|
| 683 | {
|
---|
| 684 | REPORT
|
---|
| 685 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val)
|
---|
| 686 | Throw(IndexException(m,n,*this,true));
|
---|
| 687 | return store[m*ncols_val+n];
|
---|
| 688 | }
|
---|
| 689 |
|
---|
| 690 | Real Matrix::element(int m, int n) const
|
---|
| 691 | {
|
---|
| 692 | REPORT
|
---|
| 693 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val)
|
---|
| 694 | Throw(IndexException(m,n,*this,true));
|
---|
| 695 | return store[m*ncols_val+n];
|
---|
| 696 | }
|
---|
| 697 |
|
---|
| 698 | Real& SymmetricMatrix::element(int m, int n)
|
---|
| 699 | {
|
---|
| 700 | REPORT
|
---|
| 701 | if (m<0 || n<0 || m >= nrows_val || n>=ncols_val)
|
---|
| 702 | Throw(IndexException(m,n,*this,true));
|
---|
| 703 | if (m>=n) return store[tristore(m)+n];
|
---|
| 704 | else return store[tristore(n)+m];
|
---|
| 705 | }
|
---|
| 706 |
|
---|
| 707 | Real SymmetricMatrix::element(int m, int n) const
|
---|
| 708 | {
|
---|
| 709 | REPORT
|
---|
| 710 | if (m<0 || n<0 || m >= nrows_val || n>=ncols_val)
|
---|
| 711 | Throw(IndexException(m,n,*this,true));
|
---|
| 712 | if (m>=n) return store[tristore(m)+n];
|
---|
| 713 | else return store[tristore(n)+m];
|
---|
| 714 | }
|
---|
| 715 |
|
---|
| 716 | Real& UpperTriangularMatrix::element(int m, int n)
|
---|
| 717 | {
|
---|
| 718 | REPORT
|
---|
| 719 | if (m<0 || n<m || n>=ncols_val)
|
---|
| 720 | Throw(IndexException(m,n,*this,true));
|
---|
| 721 | return store[m*ncols_val+n-tristore(m)];
|
---|
| 722 | }
|
---|
| 723 |
|
---|
| 724 | Real UpperTriangularMatrix::element(int m, int n) const
|
---|
| 725 | {
|
---|
| 726 | REPORT
|
---|
| 727 | if (m<0 || n<m || n>=ncols_val)
|
---|
| 728 | Throw(IndexException(m,n,*this,true));
|
---|
| 729 | return store[m*ncols_val+n-tristore(m)];
|
---|
| 730 | }
|
---|
| 731 |
|
---|
| 732 | Real& LowerTriangularMatrix::element(int m, int n)
|
---|
| 733 | {
|
---|
| 734 | REPORT
|
---|
| 735 | if (n<0 || m<n || m>=nrows_val)
|
---|
| 736 | Throw(IndexException(m,n,*this,true));
|
---|
| 737 | return store[tristore(m)+n];
|
---|
| 738 | }
|
---|
| 739 |
|
---|
| 740 | Real LowerTriangularMatrix::element(int m, int n) const
|
---|
| 741 | {
|
---|
| 742 | REPORT
|
---|
| 743 | if (n<0 || m<n || m>=nrows_val)
|
---|
| 744 | Throw(IndexException(m,n,*this,true));
|
---|
| 745 | return store[tristore(m)+n];
|
---|
| 746 | }
|
---|
| 747 |
|
---|
| 748 | Real& DiagonalMatrix::element(int m, int n)
|
---|
| 749 | {
|
---|
| 750 | REPORT
|
---|
| 751 | if (n<0 || m!=n || m>=nrows_val || n>=ncols_val)
|
---|
| 752 | Throw(IndexException(m,n,*this,true));
|
---|
| 753 | return store[n];
|
---|
| 754 | }
|
---|
| 755 |
|
---|
| 756 | Real DiagonalMatrix::element(int m, int n) const
|
---|
| 757 | {
|
---|
| 758 | REPORT
|
---|
| 759 | if (n<0 || m!=n || m>=nrows_val || n>=ncols_val)
|
---|
| 760 | Throw(IndexException(m,n,*this,true));
|
---|
| 761 | return store[n];
|
---|
| 762 | }
|
---|
| 763 |
|
---|
| 764 | Real& DiagonalMatrix::element(int m)
|
---|
| 765 | {
|
---|
| 766 | REPORT
|
---|
| 767 | if (m<0 || m>=nrows_val) Throw(IndexException(m,*this,true));
|
---|
| 768 | return store[m];
|
---|
| 769 | }
|
---|
| 770 |
|
---|
| 771 | Real DiagonalMatrix::element(int m) const
|
---|
| 772 | {
|
---|
| 773 | REPORT
|
---|
| 774 | if (m<0 || m>=nrows_val) Throw(IndexException(m,*this,true));
|
---|
| 775 | return store[m];
|
---|
| 776 | }
|
---|
| 777 |
|
---|
| 778 | Real& ColumnVector::element(int m)
|
---|
| 779 | {
|
---|
| 780 | REPORT
|
---|
| 781 | if (m<0 || m>= nrows_val) Throw(IndexException(m,*this,true));
|
---|
| 782 | return store[m];
|
---|
| 783 | }
|
---|
| 784 |
|
---|
| 785 | Real ColumnVector::element(int m) const
|
---|
| 786 | {
|
---|
| 787 | REPORT
|
---|
| 788 | if (m<0 || m>= nrows_val) Throw(IndexException(m,*this,true));
|
---|
| 789 | return store[m];
|
---|
| 790 | }
|
---|
| 791 |
|
---|
| 792 | Real& RowVector::element(int n)
|
---|
| 793 | {
|
---|
| 794 | REPORT
|
---|
| 795 | if (n<0 || n>= ncols_val) Throw(IndexException(n,*this,true));
|
---|
| 796 | return store[n];
|
---|
| 797 | }
|
---|
| 798 |
|
---|
| 799 | Real RowVector::element(int n) const
|
---|
| 800 | {
|
---|
| 801 | REPORT
|
---|
| 802 | if (n<0 || n>= ncols_val) Throw(IndexException(n,*this,true));
|
---|
| 803 | return store[n];
|
---|
| 804 | }
|
---|
| 805 |
|
---|
| 806 | Real& BandMatrix::element(int m, int n)
|
---|
| 807 | {
|
---|
| 808 | REPORT
|
---|
| 809 | int w = upper_val+lower_val+1; int i = lower_val+n-m;
|
---|
| 810 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
|
---|
| 811 | Throw(IndexException(m,n,*this,true));
|
---|
| 812 | return store[w*m+i];
|
---|
| 813 | }
|
---|
| 814 |
|
---|
| 815 | Real BandMatrix::element(int m, int n) const
|
---|
| 816 | {
|
---|
| 817 | REPORT
|
---|
| 818 | int w = upper_val+lower_val+1; int i = lower_val+n-m;
|
---|
| 819 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
|
---|
| 820 | Throw(IndexException(m,n,*this,true));
|
---|
| 821 | return store[w*m+i];
|
---|
| 822 | }
|
---|
| 823 |
|
---|
| 824 | Real& UpperBandMatrix::element(int m, int n)
|
---|
| 825 | {
|
---|
| 826 | REPORT
|
---|
| 827 | int w = upper_val+1; int i = n-m;
|
---|
| 828 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
|
---|
| 829 | Throw(IndexException(m,n,*this,true));
|
---|
| 830 | return store[w*m+i];
|
---|
| 831 | }
|
---|
| 832 |
|
---|
| 833 | Real UpperBandMatrix::element(int m, int n) const
|
---|
| 834 | {
|
---|
| 835 | REPORT
|
---|
| 836 | int w = upper_val+1; int i = n-m;
|
---|
| 837 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
|
---|
| 838 | Throw(IndexException(m,n,*this,true));
|
---|
| 839 | return store[w*m+i];
|
---|
| 840 | }
|
---|
| 841 |
|
---|
| 842 | Real& LowerBandMatrix::element(int m, int n)
|
---|
| 843 | {
|
---|
| 844 | REPORT
|
---|
| 845 | int w = lower_val+1; int i = lower_val+n-m;
|
---|
| 846 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
|
---|
| 847 | Throw(IndexException(m,n,*this,true));
|
---|
| 848 | return store[w*m+i];
|
---|
| 849 | }
|
---|
| 850 |
|
---|
| 851 | Real LowerBandMatrix::element(int m, int n) const
|
---|
| 852 | {
|
---|
| 853 | REPORT
|
---|
| 854 | int w = lower_val+1; int i = lower_val+n-m;
|
---|
| 855 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
|
---|
| 856 | Throw(IndexException(m,n,*this,true));
|
---|
| 857 | return store[w*m+i];
|
---|
| 858 | }
|
---|
| 859 |
|
---|
| 860 | Real& SymmetricBandMatrix::element(int m, int n)
|
---|
| 861 | {
|
---|
| 862 | REPORT
|
---|
| 863 | int w = lower_val+1;
|
---|
| 864 | if (m>=n)
|
---|
| 865 | {
|
---|
| 866 | REPORT
|
---|
| 867 | int i = lower_val+n-m;
|
---|
| 868 | if ( m>=nrows_val || n<0 || i<0 )
|
---|
| 869 | Throw(IndexException(m,n,*this,true));
|
---|
| 870 | return store[w*m+i];
|
---|
| 871 | }
|
---|
| 872 | else
|
---|
| 873 | {
|
---|
| 874 | REPORT
|
---|
| 875 | int i = lower_val+m-n;
|
---|
| 876 | if ( n>=nrows_val || m<0 || i<0 )
|
---|
| 877 | Throw(IndexException(m,n,*this,true));
|
---|
| 878 | return store[w*n+i];
|
---|
| 879 | }
|
---|
| 880 | }
|
---|
| 881 |
|
---|
| 882 | Real SymmetricBandMatrix::element(int m, int n) const
|
---|
| 883 | {
|
---|
| 884 | REPORT
|
---|
| 885 | int w = lower_val+1;
|
---|
| 886 | if (m>=n)
|
---|
| 887 | {
|
---|
| 888 | REPORT
|
---|
| 889 | int i = lower_val+n-m;
|
---|
| 890 | if ( m>=nrows_val || n<0 || i<0 )
|
---|
| 891 | Throw(IndexException(m,n,*this,true));
|
---|
| 892 | return store[w*m+i];
|
---|
| 893 | }
|
---|
| 894 | else
|
---|
| 895 | {
|
---|
| 896 | REPORT
|
---|
| 897 | int i = lower_val+m-n;
|
---|
| 898 | if ( n>=nrows_val || m<0 || i<0 )
|
---|
| 899 | Throw(IndexException(m,n,*this,true));
|
---|
| 900 | return store[w*n+i];
|
---|
| 901 | }
|
---|
| 902 | }
|
---|
| 903 |
|
---|
| 904 | #ifdef use_namespace
|
---|
| 905 | }
|
---|
| 906 | #endif
|
---|
| 907 |
|
---|
| 908 |
|
---|
| 909 | ///}
|
---|