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