newmat3.cpp
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "include.h"
00008
00009 #include "newmat.h"
00010 #include "newmatrc.h"
00011
00012 #ifdef use_namespace
00013 namespace NEWMAT {
00014 #endif
00015
00016
00017 #ifdef DO_REPORT
00018 #define REPORT { static ExeCounter ExeCount(__LINE__,3); ++ExeCount; }
00019 #else
00020 #define REPORT {}
00021 #endif
00022
00023
00024
00025
00026 #define MONITOR(what,store,storage) {}
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 void GeneralMatrix::NextRow(MatrixRowCol& mrc)
00082 {
00083 REPORT
00084 if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreRow(mrc); }
00085 mrc.rowcol++;
00086 if (mrc.rowcol<nrows_value) { REPORT this->GetRow(mrc); }
00087 else { REPORT mrc.cw -= StoreOnExit; }
00088 }
00089
00090 void GeneralMatrix::NextCol(MatrixRowCol& mrc)
00091 {
00092 REPORT
00093 if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
00094 mrc.rowcol++;
00095 if (mrc.rowcol<ncols_value) { REPORT this->GetCol(mrc); }
00096 else { REPORT mrc.cw -= StoreOnExit; }
00097 }
00098
00099 void GeneralMatrix::NextCol(MatrixColX& mrc)
00100 {
00101 REPORT
00102 if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
00103 mrc.rowcol++;
00104 if (mrc.rowcol<ncols_value) { REPORT this->GetCol(mrc); }
00105 else { REPORT mrc.cw -= StoreOnExit; }
00106 }
00107
00108
00109
00110
00111 void Matrix::GetRow(MatrixRowCol& mrc)
00112 {
00113 REPORT
00114 mrc.skip=0; mrc.storage=mrc.length=ncols_value;
00115 mrc.data=store+mrc.rowcol*ncols_value;
00116 }
00117
00118
00119 void Matrix::GetCol(MatrixRowCol& mrc)
00120 {
00121 REPORT
00122 mrc.skip=0; mrc.storage=mrc.length=nrows_value;
00123 if ( ncols_value==1 && !(mrc.cw*StoreHere) )
00124 { REPORT mrc.data=store; }
00125 else
00126 {
00127 Real* ColCopy;
00128 if ( !(mrc.cw*(HaveStore+StoreHere)) )
00129 {
00130 REPORT
00131 ColCopy = new Real [nrows_value]; MatrixErrorNoSpace(ColCopy);
00132 MONITOR_REAL_NEW("Make (MatGetCol)",nrows_value,ColCopy)
00133 mrc.data = ColCopy; mrc.cw += HaveStore;
00134 }
00135 else { REPORT ColCopy = mrc.data; }
00136 if (+(mrc.cw*LoadOnEntry))
00137 {
00138 REPORT
00139 Real* Mstore = store+mrc.rowcol; int i=nrows_value;
00140
00141 if (i) for (;;)
00142 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_value; }
00143 }
00144 }
00145 }
00146
00147 void Matrix::GetCol(MatrixColX& mrc)
00148 {
00149 REPORT
00150 mrc.skip=0; mrc.storage=nrows_value; mrc.length=nrows_value;
00151 if (+(mrc.cw*LoadOnEntry))
00152 {
00153 REPORT Real* ColCopy = mrc.data;
00154 Real* Mstore = store+mrc.rowcol; int i=nrows_value;
00155
00156 if (i) for (;;)
00157 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_value; }
00158 }
00159 }
00160
00161 void Matrix::RestoreCol(MatrixRowCol& mrc)
00162 {
00163
00164 REPORT
00165 if (+(mrc.cw*HaveStore))
00166 {
00167 REPORT
00168 Real* Mstore = store+mrc.rowcol; int i=nrows_value;
00169 Real* Cstore = mrc.data;
00170
00171 if (i) for (;;)
00172 { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_value; }
00173 }
00174 }
00175
00176 void Matrix::RestoreCol(MatrixColX& mrc)
00177 {
00178 REPORT
00179 Real* Mstore = store+mrc.rowcol; int i=nrows_value; Real* Cstore = mrc.data;
00180
00181 if (i) for (;;)
00182 { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_value; }
00183 }
00184
00185 void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }
00186
00187 void Matrix::NextCol(MatrixRowCol& mrc)
00188 {
00189 REPORT
00190 if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
00191 mrc.rowcol++;
00192 if (mrc.rowcol<ncols_value)
00193 {
00194 if (+(mrc.cw*LoadOnEntry))
00195 {
00196 REPORT
00197 Real* ColCopy = mrc.data;
00198 Real* Mstore = store+mrc.rowcol; int i=nrows_value;
00199
00200 if (i) for (;;)
00201 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_value; }
00202 }
00203 }
00204 else { REPORT mrc.cw -= StoreOnExit; }
00205 }
00206
00207 void Matrix::NextCol(MatrixColX& mrc)
00208 {
00209 REPORT
00210 if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
00211 mrc.rowcol++;
00212 if (mrc.rowcol<ncols_value)
00213 {
00214 if (+(mrc.cw*LoadOnEntry))
00215 {
00216 REPORT
00217 Real* ColCopy = mrc.data;
00218 Real* Mstore = store+mrc.rowcol; int i=nrows_value;
00219
00220 if (i) for (;;)
00221 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_value; }
00222 }
00223 }
00224 else { REPORT mrc.cw -= StoreOnExit; }
00225 }
00226
00227
00228
00229 void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
00230 {
00231 REPORT
00232 mrc.skip=mrc.rowcol; mrc.storage=1;
00233 mrc.data=store+mrc.skip; mrc.length=ncols_value;
00234 }
00235
00236 void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
00237 {
00238 REPORT
00239 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_value;
00240 if (+(mrc.cw*StoreHere))
00241 Throw(InternalException("DiagonalMatrix::GetCol(MatrixRowCol&)"));
00242 else { REPORT mrc.data=store+mrc.skip; }
00243
00244 }
00245
00246 void DiagonalMatrix::GetCol(MatrixColX& mrc)
00247 {
00248 REPORT
00249 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_value;
00250 mrc.data = mrc.store+mrc.skip;
00251 *(mrc.data)=*(store+mrc.skip);
00252 }
00253
00254 void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
00255
00256
00257 void DiagonalMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
00258
00259
00260 void DiagonalMatrix::NextCol(MatrixColX& mrc)
00261 {
00262 REPORT
00263 if (+(mrc.cw*StoreOnExit))
00264 { REPORT *(store+mrc.rowcol)=*(mrc.data); }
00265 mrc.IncrDiag();
00266 int t1 = +(mrc.cw*LoadOnEntry);
00267 if (t1 && mrc.rowcol < ncols_value)
00268 { REPORT *(mrc.data)=*(store+mrc.rowcol); }
00269 }
00270
00271
00272
00273 void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
00274 {
00275 REPORT
00276 int row = mrc.rowcol; mrc.skip=row; mrc.length=ncols_value;
00277 mrc.storage=ncols_value-row; mrc.data=store+(row*(2*ncols_value-row+1))/2;
00278 }
00279
00280
00281 void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
00282 {
00283 REPORT
00284 mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
00285 mrc.length=nrows_value; Real* ColCopy;
00286 if ( !(mrc.cw*(StoreHere+HaveStore)) )
00287 {
00288 REPORT
00289 ColCopy = new Real [nrows_value]; MatrixErrorNoSpace(ColCopy);
00290 MONITOR_REAL_NEW("Make (UT GetCol)",nrows_value,ColCopy)
00291 mrc.data = ColCopy; mrc.cw += HaveStore;
00292 }
00293 else { REPORT ColCopy = mrc.data; }
00294 if (+(mrc.cw*LoadOnEntry))
00295 {
00296 REPORT
00297 Real* Mstore = store+mrc.rowcol; int j = ncols_value;
00298
00299 if (i) for (;;)
00300 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += --j; }
00301 }
00302 }
00303
00304 void UpperTriangularMatrix::GetCol(MatrixColX& mrc)
00305 {
00306 REPORT
00307 mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
00308 mrc.length=nrows_value;
00309 if (+(mrc.cw*LoadOnEntry))
00310 {
00311 REPORT
00312 Real* ColCopy = mrc.data;
00313 Real* Mstore = store+mrc.rowcol; int j = ncols_value;
00314
00315 if (i) for (;;)
00316 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += --j; }
00317 }
00318 }
00319
00320 void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
00321 {
00322 REPORT
00323 Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols_value;
00324 Real* Cstore = mrc.data;
00325
00326 if (i) for (;;)
00327 { *Mstore = *Cstore++; if (!(--i)) break; Mstore += --j; }
00328 }
00329
00330 void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
00331
00332
00333
00334
00335 void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
00336 {
00337 REPORT
00338 int row=mrc.rowcol; mrc.skip=0; mrc.storage=row+1; mrc.length=ncols_value;
00339 mrc.data=store+(row*(row+1))/2;
00340 }
00341
00342 void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
00343 {
00344 REPORT
00345 int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows_value;
00346 int i=nrows_value-col; mrc.storage=i; Real* ColCopy;
00347 if ( +(mrc.cw*(StoreHere+HaveStore)) )
00348 { REPORT ColCopy = mrc.data; }
00349 else
00350 {
00351 REPORT
00352 ColCopy = new Real [nrows_value]; MatrixErrorNoSpace(ColCopy);
00353 MONITOR_REAL_NEW("Make (LT GetCol)",nrows_value,ColCopy)
00354 mrc.cw += HaveStore; mrc.data = ColCopy;
00355 }
00356
00357 if (+(mrc.cw*LoadOnEntry))
00358 {
00359 REPORT
00360 Real* Mstore = store+(col*(col+3))/2;
00361
00362 if (i) for (;;)
00363 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00364 }
00365 }
00366
00367 void LowerTriangularMatrix::GetCol(MatrixColX& mrc)
00368 {
00369 REPORT
00370 int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows_value;
00371 int i=nrows_value-col; mrc.storage=i; mrc.data = mrc.store + col;
00372
00373 if (+(mrc.cw*LoadOnEntry))
00374 {
00375 REPORT Real* ColCopy = mrc.data;
00376 Real* Mstore = store+(col*(col+3))/2;
00377
00378 if (i) for (;;)
00379 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00380 }
00381 }
00382
00383 void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
00384 {
00385 REPORT
00386 int col=mrc.rowcol; Real* Cstore = mrc.data;
00387 Real* Mstore = store+(col*(col+3))/2; int i=nrows_value-col;
00388
00389 if (i) for (;;)
00390 { *Mstore = *Cstore++; if (!(--i)) break; Mstore += ++col; }
00391 }
00392
00393 void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
00394
00395
00396
00397
00398 void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
00399 {
00400 REPORT
00401 mrc.skip=0; int row=mrc.rowcol; mrc.length=ncols_value;
00402 if (+(mrc.cw*DirectPart))
00403 { REPORT mrc.storage=row+1; mrc.data=store+(row*(row+1))/2; }
00404 else
00405 {
00406
00407 if (+(mrc.cw*StoreOnExit))
00408 Throw(InternalException("SymmetricMatrix::GetRow(MatrixRowCol&)"));
00409 mrc.storage=ncols_value; Real* RowCopy;
00410 if (!(mrc.cw*HaveStore))
00411 {
00412 REPORT
00413 RowCopy = new Real [ncols_value]; MatrixErrorNoSpace(RowCopy);
00414 MONITOR_REAL_NEW("Make (SymGetRow)",ncols_value,RowCopy)
00415 mrc.cw += HaveStore; mrc.data = RowCopy;
00416 }
00417 else { REPORT RowCopy = mrc.data; }
00418 if (+(mrc.cw*LoadOnEntry))
00419 {
00420 REPORT
00421 Real* Mstore = store+(row*(row+1))/2; int i = row;
00422 while (i--) *RowCopy++ = *Mstore++;
00423 i = ncols_value-row;
00424
00425 if (i) for (;;)
00426 { *RowCopy++ = *Mstore; if (!(--i)) break; Mstore += ++row; }
00427 }
00428 }
00429 }
00430
00431 void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
00432 {
00433
00434 if (+(mrc.cw*StoreHere))
00435 Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
00436
00437 int col=mrc.rowcol; mrc.length=nrows_value;
00438 REPORT
00439 mrc.skip=0;
00440 if (+(mrc.cw*DirectPart))
00441 { REPORT mrc.storage=col+1; mrc.data=store+(col*(col+1))/2; }
00442 else
00443 {
00444
00445 if (+(mrc.cw*StoreOnExit))
00446 Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
00447
00448 mrc.storage=ncols_value; Real* ColCopy;
00449 if ( +(mrc.cw*HaveStore)) { REPORT ColCopy = mrc.data; }
00450 else
00451 {
00452 REPORT
00453 ColCopy = new Real [ncols_value]; MatrixErrorNoSpace(ColCopy);
00454 MONITOR_REAL_NEW("Make (SymGetCol)",ncols_value,ColCopy)
00455 mrc.cw += HaveStore; mrc.data = ColCopy;
00456 }
00457 if (+(mrc.cw*LoadOnEntry))
00458 {
00459 REPORT
00460 Real* Mstore = store+(col*(col+1))/2; int i = col;
00461 while (i--) *ColCopy++ = *Mstore++;
00462 i = ncols_value-col;
00463
00464 if (i) for (;;)
00465 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00466 }
00467 }
00468 }
00469
00470 void SymmetricMatrix::GetCol(MatrixColX& mrc)
00471 {
00472 int col=mrc.rowcol; mrc.length=nrows_value;
00473 if (+(mrc.cw*DirectPart))
00474 {
00475 REPORT
00476 mrc.skip=col; int i=nrows_value-col; mrc.storage=i;
00477 mrc.data = mrc.store+col;
00478 if (+(mrc.cw*LoadOnEntry))
00479 {
00480 REPORT
00481 Real* ColCopy = mrc.data;
00482 Real* Mstore = store+(col*(col+3))/2;
00483
00484 if (i) for (;;)
00485 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00486 }
00487 }
00488 else
00489 {
00490 REPORT
00491
00492 if (+(mrc.cw*StoreOnExit))
00493 Throw(InternalException("SymmetricMatrix::GetCol(MatrixColX&)"));
00494
00495 mrc.skip=0; mrc.storage=ncols_value;
00496 if (+(mrc.cw*LoadOnEntry))
00497 {
00498 REPORT
00499 Real* ColCopy = mrc.data;
00500 Real* Mstore = store+(col*(col+1))/2; int i = col;
00501 while (i--) *ColCopy++ = *Mstore++;
00502 i = ncols_value-col;
00503
00504 if (i) for (;;)
00505 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00506 }
00507 }
00508 }
00509
00510
00511
00512 void SymmetricMatrix::RestoreCol(MatrixColX& mrc)
00513 {
00514 REPORT
00515
00516 int col=mrc.rowcol; Real* Cstore = mrc.data;
00517 Real* Mstore = store+(col*(col+3))/2; int i = nrows_value-col;
00518
00519 if (i) for (;;)
00520 { *Mstore = *Cstore++; if (!(--i)) break; Mstore+= ++col; }
00521
00522 }
00523
00524
00525
00526 void RowVector::GetCol(MatrixRowCol& mrc)
00527 {
00528 REPORT
00529
00530 if (+(mrc.cw*StoreHere))
00531 Throw(InternalException("RowVector::GetCol(MatrixRowCol&)"));
00532
00533 mrc.skip=0; mrc.storage=1; mrc.length=nrows_value;
00534 mrc.data = store+mrc.rowcol;
00535
00536 }
00537
00538 void RowVector::GetCol(MatrixColX& mrc)
00539 {
00540 REPORT
00541 mrc.skip=0; mrc.storage=1; mrc.length=nrows_value;
00542 if (mrc.cw >= LoadOnEntry)
00543 { REPORT *(mrc.data) = *(store+mrc.rowcol); }
00544
00545 }
00546
00547 void RowVector::NextCol(MatrixRowCol& mrc)
00548 { REPORT mrc.rowcol++; mrc.data++; }
00549
00550 void RowVector::NextCol(MatrixColX& mrc)
00551 {
00552 if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.data); }
00553
00554 mrc.rowcol++;
00555 if (mrc.rowcol < ncols_value)
00556 {
00557 if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.data)=*(store+mrc.rowcol); }
00558 }
00559 else { REPORT mrc.cw -= StoreOnExit; }
00560 }
00561
00562 void RowVector::RestoreCol(MatrixColX& mrc)
00563 { REPORT *(store+mrc.rowcol)=*(mrc.data); }
00564
00565
00566
00567
00568 void BandMatrix::GetRow(MatrixRowCol& mrc)
00569 {
00570 REPORT
00571 int r = mrc.rowcol; int w = lower+1+upper; mrc.length=ncols_value;
00572 int s = r-lower;
00573 if (s<0) { mrc.data = store+(r*w-s); w += s; s = 0; }
00574 else mrc.data = store+r*w;
00575 mrc.skip = s; s += w-ncols_value; if (s>0) w -= s; mrc.storage = w;
00576 }
00577
00578
00579
00580 void BandMatrix::NextRow(MatrixRowCol& mrc)
00581 {
00582 REPORT
00583 int r = ++mrc.rowcol;
00584 if (r<=lower) { mrc.storage++; mrc.data += lower+upper; }
00585 else { mrc.skip++; mrc.data += lower+upper+1; }
00586 if (r>=ncols_value-upper) mrc.storage--;
00587 }
00588
00589 void BandMatrix::GetCol(MatrixRowCol& mrc)
00590 {
00591 REPORT
00592 int c = mrc.rowcol; int n = lower+upper; int w = n+1;
00593 mrc.length=nrows_value; Real* ColCopy;
00594 int b; int s = c-upper;
00595 if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
00596 mrc.skip = s; s += w-nrows_value; if (s>0) w -= s; mrc.storage = w;
00597 if ( +(mrc.cw*(StoreHere+HaveStore)) )
00598 { REPORT ColCopy = mrc.data; }
00599 else
00600 {
00601 REPORT
00602 ColCopy = new Real [n+1]; MatrixErrorNoSpace(ColCopy);
00603 MONITOR_REAL_NEW("Make (BMGetCol)",n+1,ColCopy)
00604 mrc.cw += HaveStore; mrc.data = ColCopy;
00605 }
00606
00607 if (+(mrc.cw*LoadOnEntry))
00608 {
00609 REPORT
00610 Real* Mstore = store+b;
00611
00612 if (w) for (;;)
00613 { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
00614 }
00615 }
00616
00617 void BandMatrix::GetCol(MatrixColX& mrc)
00618 {
00619 REPORT
00620 int c = mrc.rowcol; int n = lower+upper; int w = n+1;
00621 mrc.length=nrows_value; int b; int s = c-upper;
00622 if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
00623 mrc.skip = s; s += w-nrows_value; if (s>0) w -= s; mrc.storage = w;
00624 mrc.data = mrc.store+mrc.skip;
00625
00626 if (+(mrc.cw*LoadOnEntry))
00627 {
00628 REPORT
00629 Real* ColCopy = mrc.data; Real* Mstore = store+b;
00630
00631 if (w) for (;;)
00632 { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
00633 }
00634 }
00635
00636 void BandMatrix::RestoreCol(MatrixRowCol& mrc)
00637 {
00638 REPORT
00639 int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
00640 Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
00641 Real* Cstore = mrc.data;
00642 int w = mrc.storage;
00643
00644 if (w) for (;;)
00645 { *Mstore = *Cstore++; if (!(--w)) break; Mstore += n; }
00646 }
00647
00648
00649
00650 void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
00651 {
00652 REPORT
00653 int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
00654 mrc.length = ncols_value;
00655 if (s<0) { w1 += s; o -= s; s = 0; }
00656 mrc.skip = s;
00657
00658 if (+(mrc.cw*DirectPart))
00659 { REPORT mrc.data = store+o; mrc.storage = w1; }
00660 else
00661 {
00662
00663 if (+(mrc.cw*StoreOnExit))
00664 Throw(InternalException("SymmetricBandMatrix::GetRow(MatrixRowCol&)"));
00665 int w = w1+lower; s += w-ncols_value; Real* RowCopy;
00666 if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
00667 if (!(mrc.cw*HaveStore))
00668 {
00669 REPORT
00670 RowCopy = new Real [2*lower+1]; MatrixErrorNoSpace(RowCopy);
00671 MONITOR_REAL_NEW("Make (SmBGetRow)",2*lower+1,RowCopy)
00672 mrc.cw += HaveStore; mrc.data = RowCopy;
00673 }
00674 else { REPORT RowCopy = mrc.data; }
00675
00676 if (+(mrc.cw*LoadOnEntry) && ncols_value > 0)
00677 {
00678 REPORT
00679 Real* Mstore = store+o;
00680 while (w1--) *RowCopy++ = *Mstore++;
00681 Mstore--;
00682 while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
00683 }
00684 }
00685 }
00686
00687 void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
00688 {
00689
00690 if (+(mrc.cw*StoreHere))
00691 Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
00692
00693 int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows_value;
00694 REPORT
00695 int s = c-lower; int o = c*w1;
00696 if (s<0) { w1 += s; o -= s; s = 0; }
00697 mrc.skip = s;
00698
00699 if (+(mrc.cw*DirectPart))
00700 { REPORT mrc.data = store+o; mrc.storage = w1; }
00701 else
00702 {
00703
00704 if (+(mrc.cw*StoreOnExit))
00705 Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
00706 int w = w1+lower; s += w-ncols_value; Real* ColCopy;
00707 if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
00708
00709 if ( +(mrc.cw*HaveStore) ) { REPORT ColCopy = mrc.data; }
00710 else
00711 {
00712 REPORT ColCopy = new Real [2*lower+1]; MatrixErrorNoSpace(ColCopy);
00713 MONITOR_REAL_NEW("Make (SmBGetCol)",2*lower+1,ColCopy)
00714 mrc.cw += HaveStore; mrc.data = ColCopy;
00715 }
00716
00717 if (+(mrc.cw*LoadOnEntry))
00718 {
00719 REPORT
00720 Real* Mstore = store+o;
00721 while (w1--) *ColCopy++ = *Mstore++;
00722 Mstore--;
00723 while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
00724 }
00725 }
00726 }
00727
00728 void SymmetricBandMatrix::GetCol(MatrixColX& mrc)
00729 {
00730 int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows_value;
00731 if (+(mrc.cw*DirectPart))
00732 {
00733 REPORT
00734 int b = c*w1+lower;
00735 mrc.skip = c; c += w1-nrows_value; w1 -= c; mrc.storage = w1;
00736 Real* ColCopy = mrc.data = mrc.store+mrc.skip;
00737 if (+(mrc.cw*LoadOnEntry))
00738 {
00739 REPORT
00740 Real* Mstore = store+b;
00741
00742 if (w1) for (;;)
00743 { *ColCopy++ = *Mstore; if (!(--w1)) break; Mstore += lower; }
00744 }
00745 }
00746 else
00747 {
00748 REPORT
00749
00750 if (+(mrc.cw*StoreOnExit))
00751 Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixColX&)"));
00752 int s = c-lower; int o = c*w1;
00753 if (s<0) { w1 += s; o -= s; s = 0; }
00754 mrc.skip = s;
00755
00756 int w = w1+lower; s += w-ncols_value;
00757 if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
00758
00759 Real* ColCopy = mrc.data = mrc.store+mrc.skip;
00760
00761 if (+(mrc.cw*LoadOnEntry))
00762 {
00763 REPORT
00764 Real* Mstore = store+o;
00765 while (w1--) *ColCopy++ = *Mstore++;
00766 Mstore--;
00767 while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
00768 }
00769
00770 }
00771 }
00772
00773 void SymmetricBandMatrix::RestoreCol(MatrixColX& mrc)
00774 {
00775 REPORT
00776 int c = mrc.rowcol;
00777 Real* Mstore = store + c*lower+c+lower;
00778 Real* Cstore = mrc.data; int w = mrc.storage;
00779
00780 if (w) for (;;)
00781 { *Mstore = *Cstore++; if (!(--w)) break; Mstore += lower; }
00782 }
00783
00784
00785
00786 void IdentityMatrix::GetRow(MatrixRowCol& mrc)
00787 {
00788 REPORT
00789 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.data=store; mrc.length=ncols_value;
00790 }
00791
00792 void IdentityMatrix::GetCol(MatrixRowCol& mrc)
00793 {
00794 REPORT
00795 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_value;
00796 if (+(mrc.cw*StoreHere))
00797 Throw(InternalException("IdentityMatrix::GetCol(MatrixRowCol&)"));
00798 else { REPORT mrc.data=store; }
00799 }
00800
00801 void IdentityMatrix::GetCol(MatrixColX& mrc)
00802 {
00803 REPORT
00804 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_value;
00805 mrc.data = mrc.store+mrc.skip; *(mrc.data)=*store;
00806 }
00807
00808 void IdentityMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrId(); }
00809
00810 void IdentityMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrId(); }
00811
00812 void IdentityMatrix::NextCol(MatrixColX& mrc)
00813 {
00814 REPORT
00815 if (+(mrc.cw*StoreOnExit)) { REPORT *store=*(mrc.data); }
00816 mrc.IncrDiag();
00817 int t1 = +(mrc.cw*LoadOnEntry);
00818 if (t1 && mrc.rowcol < ncols_value) { REPORT *(mrc.data)=*store; }
00819 }
00820
00821
00822
00823
00824
00825
00826 MatrixRowCol::~MatrixRowCol()
00827 {
00828 if (+(cw*HaveStore))
00829 {
00830 MONITOR_REAL_DELETE("Free (RowCol)",-1,data)
00831 delete [] data;
00832 }
00833 }
00834
00835 MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }
00836
00837 MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
00838
00839 MatrixColX::~MatrixColX() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
00840
00841 #ifdef use_namespace
00842 }
00843 #endif
00844