NEWMAT3.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
NEWMAT3.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ newmat3.cpp        Matrix get and restore rows and columns
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #include "include.h"
00007 
00008 #include "newmat.h"
00009 #include "newmatrc.h"
00010 
00011 #ifdef use_namespace
00012 namespace NEWMAT {
00013 #endif
00014 
00015 #ifdef DO_REPORT
00016 #define REPORT { static ExeCounter ExeCount(__LINE__,3); ++ExeCount; }
00017 #else
00018 #define REPORT {}
00019 #endif
00020 
00021     //#define MONITOR(what,storage,store) { cout << what << " " << storage << " at " << (long)store << "\n"; }
00022 
00023 #define MONITOR(what,store,storage) {}
00024 
00025     // Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
00026     //
00027     // LoadOnEntry:
00028     //    Load data into MatrixRow or Col dummy array under GetRow or GetCol
00029     // StoreOnExit:
00030     //    Restore data to original matrix under RestoreRow or RestoreCol
00031     // DirectPart:
00032     //    Load or restore only part directly stored; must be set with StoreOnExit
00033     //    Still have decide how to handle this with symmetric
00034     // StoreHere:
00035     //    used in columns only - store data at supplied storage address;
00036     //    used for GetCol, NextCol & RestoreCol. No need to fill out zeros
00037     // HaveStore:
00038     //    dummy array has been assigned (internal use only).
00039 
00040     // For symmetric matrices, treat columns as rows unless StoreHere is set;
00041     // then stick to columns as this will give better performance for doing
00042     // inverses
00043 
00044     // How components are used:
00045 
00046     // Use rows wherever possible in preference to columns
00047 
00048     // Columns without StoreHere are used in in-exact transpose, sum column,
00049     // multiply a column vector, and maybe in future access to column,
00050     // additional multiply functions, add transpose
00051 
00052     // Columns with StoreHere are used in exact transpose (not symmetric matrices
00053     // or vectors, load only)
00054 
00055     // Columns with MatrixColX (Store to full column) are used in inverse and solve
00056 
00057     // Functions required for each matrix class
00058 
00059     // GetRow(MatrixRowCol& mrc)
00060     // GetCol(MatrixRowCol& mrc)
00061     // GetCol(MatrixColX& mrc)
00062     // RestoreRow(MatrixRowCol& mrc)
00063     // RestoreCol(MatrixRowCol& mrc)
00064     // RestoreCol(MatrixColX& mrc)
00065     // NextRow(MatrixRowCol& mrc)
00066     // NextCol(MatrixRowCol& mrc)
00067     // NextCol(MatrixColX& mrc)
00068 
00069     // The Restore routines assume StoreOnExit has already been checked
00070     // Defaults for the Next routines are given below
00071     // Assume can't have both !DirectPart && StoreHere for MatrixRowCol routines
00072 
00073     // Default NextRow and NextCol:
00074     // will work as a default but need to override NextRow for efficiency
00075 
00076     void GeneralMatrix::NextRow(MatrixRowCol& mrc) {
00077         REPORT
00078             if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreRow(mrc); }
00079         mrc.rowcol++;
00080         if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
00081         else { REPORT mrc.cw -= StoreOnExit; }
00082     }
00083 
00084     void GeneralMatrix::NextCol(MatrixRowCol& mrc) {
00085         REPORT                                        // 423
00086             if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
00087         mrc.rowcol++;
00088         if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
00089         else { REPORT mrc.cw -= StoreOnExit; }
00090     }
00091 
00092     void GeneralMatrix::NextCol(MatrixColX& mrc) {
00093         REPORT                                        // 423
00094             if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
00095         mrc.rowcol++;
00096         if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
00097         else { REPORT mrc.cw -= StoreOnExit; }
00098     }
00099 
00100     // routines for matrix
00101 
00102     void Matrix::GetRow(MatrixRowCol& mrc) {
00103         REPORT
00104             mrc.skip=0; mrc.storage=mrc.length=ncols; mrc.data=store+mrc.rowcol*ncols;
00105     }
00106 
00107     void Matrix::GetCol(MatrixRowCol& mrc) {
00108         REPORT
00109             mrc.skip=0; mrc.storage=mrc.length=nrows;
00110         if ( ncols==1 && !(mrc.cw*StoreHere) )        // ColumnVector
00111         { REPORT mrc.data=store; }
00112         else {
00113             Real* ColCopy;
00114             if ( !(mrc.cw*(HaveStore+StoreHere)) ) {
00115                 REPORT
00116                     ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
00117                 MONITOR_REAL_NEW("Make (MatGetCol)",nrows,ColCopy)
00118                     mrc.data = ColCopy; mrc.cw += HaveStore;
00119             }
00120             else { REPORT ColCopy = mrc.data; }
00121             if (+(mrc.cw*LoadOnEntry)) {
00122                 REPORT
00123                     Real* Mstore = store+mrc.rowcol; int i=nrows;
00124                 while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
00125             }
00126         }
00127     }
00128 
00129     void Matrix::GetCol(MatrixColX& mrc) {
00130         REPORT
00131             mrc.skip=0; mrc.storage=nrows; mrc.length=nrows;
00132             if (+(mrc.cw*LoadOnEntry)) {
00133                 REPORT  Real* ColCopy = mrc.data;
00134                 Real* Mstore = store+mrc.rowcol; int i=nrows;
00135                 while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
00136             }
00137     }
00138 
00139     void Matrix::RestoreCol(MatrixRowCol& mrc) {
00140         // always check StoreOnExit before calling RestoreCol
00141         REPORT                                        // 429
00142             if (+(mrc.cw*HaveStore)) {
00143                 REPORT                                      // 426
00144                     Real* Mstore = store+mrc.rowcol; int i=nrows;
00145                 Real* Cstore = mrc.data;
00146                 while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
00147             }
00148     }
00149 
00150     void Matrix::RestoreCol(MatrixColX& mrc) {
00151         REPORT
00152             Real* Mstore = store+mrc.rowcol; int i=nrows; Real* Cstore = mrc.data;
00153             while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
00154     }
00155 
00156     void Matrix::NextRow(MatrixRowCol& mrc) {       // 1808
00157         REPORT mrc.IncrMat();
00158     }
00159 
00160     void Matrix::NextCol(MatrixRowCol& mrc) {
00161         REPORT                                        // 632
00162             if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
00163         mrc.rowcol++;
00164         if (mrc.rowcol<ncols) {
00165             if (+(mrc.cw*LoadOnEntry)) {
00166                 REPORT
00167                     Real* ColCopy = mrc.data;
00168                 Real* Mstore = store+mrc.rowcol; int i=nrows;
00169                 while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
00170             }
00171         }
00172         else { REPORT mrc.cw -= StoreOnExit; }
00173     }
00174 
00175     void Matrix::NextCol(MatrixColX& mrc) {
00176         REPORT
00177             if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
00178         mrc.rowcol++;
00179         if (mrc.rowcol<ncols) {
00180             if (+(mrc.cw*LoadOnEntry)) {
00181                 REPORT
00182                     Real* ColCopy = mrc.data;
00183                 Real* Mstore = store+mrc.rowcol; int i=nrows;
00184                 while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
00185             }
00186         }
00187         else { REPORT mrc.cw -= StoreOnExit; }
00188     }
00189 
00190     // routines for diagonal matrix
00191 
00192     void DiagonalMatrix::GetRow(MatrixRowCol& mrc) {
00193         REPORT
00194             mrc.skip=mrc.rowcol; mrc.storage=1;
00195         mrc.data=store+mrc.skip; mrc.length=ncols;
00196     }
00197 
00198     void DiagonalMatrix::GetCol(MatrixRowCol& mrc) {
00199         REPORT
00200             mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows;
00201             if (+(mrc.cw*StoreHere))                      // shouldn't happen
00202                 Throw(InternalException("DiagonalMatrix::GetCol(MatrixRowCol&)"));
00203             else  { REPORT mrc.data=store+mrc.skip; }
00204             // not accessed
00205     }
00206 
00207     void DiagonalMatrix::GetCol(MatrixColX& mrc) {
00208         REPORT
00209             mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows;
00210             mrc.data = mrc.store+mrc.skip;
00211             *(mrc.data)=*(store+mrc.skip);
00212     }
00213 
00214     void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
00215     // 800
00216 
00217     void DiagonalMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
00218     // not accessed
00219 
00220     void DiagonalMatrix::NextCol(MatrixColX& mrc) {
00221         REPORT
00222             if (+(mrc.cw*StoreOnExit))
00223             { REPORT *(store+mrc.rowcol)=*(mrc.data); }
00224         mrc.IncrDiag();
00225         int t1 = +(mrc.cw*LoadOnEntry);
00226         if (t1 && mrc.rowcol < ncols)
00227         { REPORT *(mrc.data)=*(store+mrc.rowcol); }
00228     }
00229 
00230     // routines for upper triangular matrix
00231 
00232     void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc) {
00233         REPORT
00234             int row = mrc.rowcol; mrc.skip=row; mrc.length=ncols;
00235             mrc.storage=ncols-row; mrc.data=store+(row*(2*ncols-row+1))/2;
00236     }
00237 
00238     void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc) {
00239         REPORT
00240             mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
00241             mrc.length=nrows; Real* ColCopy;
00242             if ( !(mrc.cw*(StoreHere+HaveStore)) ) {
00243                 REPORT                                      // not accessed
00244                     ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
00245                 MONITOR_REAL_NEW("Make (UT GetCol)",nrows,ColCopy)
00246                     mrc.data = ColCopy; mrc.cw += HaveStore;
00247             }
00248             else { REPORT ColCopy = mrc.data; }
00249             if (+(mrc.cw*LoadOnEntry)) {
00250                 REPORT
00251                     Real* Mstore = store+mrc.rowcol; int j = ncols;
00252                 while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
00253             }
00254     }
00255 
00256     void UpperTriangularMatrix::GetCol(MatrixColX& mrc) {
00257         REPORT
00258             mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
00259             mrc.length=nrows;
00260             if (+(mrc.cw*LoadOnEntry)) {
00261                 REPORT
00262                     Real* ColCopy = mrc.data;
00263                 Real* Mstore = store+mrc.rowcol; int j = ncols;
00264                 while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
00265             }
00266     }
00267 
00268     void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc) {
00269         REPORT
00270             Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
00271             Real* Cstore = mrc.data;
00272             while (i--) { *Mstore = *Cstore++; Mstore += --j; }
00273     }
00274 
00275     void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
00276     // 722
00277 
00278     // routines for lower triangular matrix
00279 
00280     void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc) {
00281         REPORT
00282             int row=mrc.rowcol; mrc.skip=0; mrc.storage=row+1; mrc.length=ncols;
00283             mrc.data=store+(row*(row+1))/2;
00284     }
00285 
00286     void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc) {
00287         REPORT
00288             int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows;
00289             int i=nrows-col; mrc.storage=i; Real* ColCopy;
00290             if ( +(mrc.cw*(StoreHere+HaveStore)) )
00291             { REPORT  ColCopy = mrc.data; }
00292             else {
00293                 REPORT                                      // not accessed
00294                     ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
00295                 MONITOR_REAL_NEW("Make (LT GetCol)",nrows,ColCopy)
00296                     mrc.cw += HaveStore; mrc.data = ColCopy;
00297             }
00298 
00299             if (+(mrc.cw*LoadOnEntry)) {
00300                 REPORT
00301                     Real* Mstore = store+(col*(col+3))/2;
00302                 while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00303             }
00304     }
00305 
00306     void LowerTriangularMatrix::GetCol(MatrixColX& mrc) {
00307         REPORT
00308             int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows;
00309             int i=nrows-col; mrc.storage=i; mrc.data = mrc.store + col;
00310 
00311             if (+(mrc.cw*LoadOnEntry)) {
00312                 REPORT  Real* ColCopy = mrc.data;
00313                 Real* Mstore = store+(col*(col+3))/2;
00314                 while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00315             }
00316     }
00317 
00318     void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc) {
00319         REPORT
00320             int col=mrc.rowcol; Real* Cstore = mrc.data;
00321         Real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
00322         while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
00323     }
00324 
00325     void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
00326     //712
00327 
00328     // routines for symmetric matrix
00329 
00330     void SymmetricMatrix::GetRow(MatrixRowCol& mrc) {
00331         REPORT                                        //571
00332             mrc.skip=0; int row=mrc.rowcol; mrc.length=ncols;
00333             if (+(mrc.cw*DirectPart))
00334             { REPORT mrc.storage=row+1; mrc.data=store+(row*(row+1))/2; }
00335             else {
00336                 // don't allow StoreOnExit and !DirectPart
00337                 if (+(mrc.cw*StoreOnExit))
00338                     Throw(InternalException("SymmetricMatrix::GetRow(MatrixRowCol&)"));
00339                 mrc.storage=ncols; Real* RowCopy;
00340                 if (!(mrc.cw*HaveStore)) {
00341                     REPORT
00342                         RowCopy = new Real [ncols]; MatrixErrorNoSpace(RowCopy);
00343                     MONITOR_REAL_NEW("Make (SymGetRow)",ncols,RowCopy)
00344                         mrc.cw += HaveStore; mrc.data = RowCopy;
00345                 }
00346                 else { REPORT RowCopy = mrc.data; }
00347                 if (+(mrc.cw*LoadOnEntry)) {
00348                     REPORT                                    // 544
00349                         Real* Mstore = store+(row*(row+1))/2; int i = row;
00350                     while (i--) *RowCopy++ = *Mstore++;
00351                     i = ncols-row;
00352                     while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
00353                 }
00354             }
00355     }
00356 
00357     void SymmetricMatrix::GetCol(MatrixRowCol& mrc) {
00358         // don't allow StoreHere
00359         if (+(mrc.cw*StoreHere))
00360             Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
00361 
00362         int col=mrc.rowcol; mrc.length=nrows;
00363         REPORT
00364             mrc.skip=0;
00365         if (+(mrc.cw*DirectPart))                     // actually get row ??
00366         { REPORT mrc.storage=col+1; mrc.data=store+(col*(col+1))/2; }
00367         else {
00368             // don't allow StoreOnExit and !DirectPart
00369             if (+(mrc.cw*StoreOnExit))
00370                 Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
00371 
00372             mrc.storage=ncols; Real* ColCopy;
00373             if ( +(mrc.cw*HaveStore)) { REPORT ColCopy = mrc.data; }
00374             else {
00375                 REPORT                                    // not accessed
00376                     ColCopy = new Real [ncols]; MatrixErrorNoSpace(ColCopy);
00377                 MONITOR_REAL_NEW("Make (SymGetCol)",ncols,ColCopy)
00378                     mrc.cw += HaveStore; mrc.data = ColCopy;
00379             }
00380             if (+(mrc.cw*LoadOnEntry)) {
00381                 REPORT
00382                     Real* Mstore = store+(col*(col+1))/2; int i = col;
00383                 while (i--) *ColCopy++ = *Mstore++;
00384                 i = ncols-col;
00385                 while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00386             }
00387         }
00388     }
00389 
00390     void SymmetricMatrix::GetCol(MatrixColX& mrc) {
00391         int col=mrc.rowcol; mrc.length=nrows;
00392         if (+(mrc.cw*DirectPart)) {
00393             REPORT
00394                 mrc.skip=col; int i=nrows-col; mrc.storage=i;
00395                 mrc.data = mrc.store+col;
00396                 if (+(mrc.cw*LoadOnEntry)) {
00397                     REPORT                                    // not accessed
00398                         Real* ColCopy = mrc.data;
00399                     Real* Mstore = store+(col*(col+3))/2;
00400                     while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00401                 }
00402         }
00403         else {
00404             REPORT
00405                 // don't allow StoreOnExit and !DirectPart
00406                 if (+(mrc.cw*StoreOnExit))
00407                     Throw(InternalException("SymmetricMatrix::GetCol(MatrixColX&)"));
00408 
00409             mrc.skip=0; mrc.storage=ncols;
00410             if (+(mrc.cw*LoadOnEntry)) {
00411                 REPORT
00412                     Real* ColCopy = mrc.data;
00413                 Real* Mstore = store+(col*(col+1))/2; int i = col;
00414                 while (i--) *ColCopy++ = *Mstore++;
00415                 i = ncols-col;
00416                 while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00417             }
00418         }
00419     }
00420 
00421     // Don't need RestoreRow because we don't allow !DirectPart && StoreOnExit
00422 
00423     void SymmetricMatrix::RestoreCol(MatrixColX& mrc) {
00424         REPORT
00425             // Really do restore column
00426             int col=mrc.rowcol; Real* Cstore = mrc.data;
00427         Real* Mstore = store+(col*(col+3))/2; int i = nrows-col;
00428         while (i--) { *Mstore = *Cstore++; Mstore+= ++col; }
00429 
00430     }
00431 
00432     // routines for row vector
00433 
00434     void RowVector::GetCol(MatrixRowCol& mrc) {
00435         REPORT
00436             // don't allow StoreHere
00437             if (+(mrc.cw*StoreHere))
00438                 Throw(InternalException("RowVector::GetCol(MatrixRowCol&)"));
00439 
00440         mrc.skip=0; mrc.storage=1; mrc.length=nrows; mrc.data = store+mrc.rowcol;
00441 
00442     }
00443 
00444     void RowVector::GetCol(MatrixColX& mrc) {
00445         REPORT
00446             mrc.skip=0; mrc.storage=1; mrc.length=nrows;
00447             if (mrc.cw >= LoadOnEntry)
00448             { REPORT *(mrc.data) = *(store+mrc.rowcol); }
00449 
00450     }
00451 
00452     void RowVector::NextCol(MatrixRowCol& mrc)
00453     { REPORT mrc.rowcol++; mrc.data++; }
00454 
00455     void RowVector::NextCol(MatrixColX& mrc) {
00456         if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.data); }
00457 
00458         mrc.rowcol++;
00459         if (mrc.rowcol < ncols) {
00460             if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.data)=*(store+mrc.rowcol); }
00461         }
00462         else { REPORT mrc.cw -= StoreOnExit; }
00463     }
00464 
00465     void RowVector::RestoreCol(MatrixColX& mrc) {
00466         // not accessed
00467         REPORT *(store+mrc.rowcol)=*(mrc.data);
00468     }
00469 
00470     // routines for band matrices
00471 
00472     void BandMatrix::GetRow(MatrixRowCol& mrc) {
00473         REPORT
00474             int r = mrc.rowcol; int w = lower+1+upper; mrc.length=ncols;
00475             int s = r-lower;
00476             if (s<0) { mrc.data = store+(r*w-s); w += s; s = 0; }
00477             else mrc.data = store+r*w;
00478             mrc.skip = s; s += w-ncols; if (s>0) w -= s; mrc.storage = w;
00479     }
00480 
00481     // should make special versions of this for upper and lower band matrices
00482 
00483     void BandMatrix::NextRow(MatrixRowCol& mrc) {
00484         REPORT
00485             int r = ++mrc.rowcol;
00486         if (r<=lower) { mrc.storage++; mrc.data += lower+upper; }
00487         else  { mrc.skip++; mrc.data += lower+upper+1; }
00488         if (r>=ncols-upper) mrc.storage--;
00489     }
00490 
00491     void BandMatrix::GetCol(MatrixRowCol& mrc) {
00492         REPORT
00493             int c = mrc.rowcol; int n = lower+upper; int w = n+1;
00494             mrc.length=nrows; Real* ColCopy;
00495             int b; int s = c-upper;
00496             if (s<=0) { w += s; s = 0; b = c+lower; }
00497             else b = s*w+n;
00498             mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
00499             if ( +(mrc.cw*(StoreHere+HaveStore)) )
00500             { REPORT ColCopy = mrc.data; }
00501             else {
00502                 REPORT
00503                     ColCopy = new Real [n+1]; MatrixErrorNoSpace(ColCopy);
00504                 MONITOR_REAL_NEW("Make (BMGetCol)",n+1,ColCopy)
00505                     mrc.cw += HaveStore; mrc.data = ColCopy;
00506             }
00507 
00508             if (+(mrc.cw*LoadOnEntry)) {
00509                 REPORT
00510                     Real* Mstore = store+b;
00511                 while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
00512             }
00513     }
00514 
00515     void BandMatrix::GetCol(MatrixColX& mrc) {
00516         REPORT
00517             int c = mrc.rowcol; int n = lower+upper; int w = n+1;
00518             mrc.length=nrows; int b; int s = c-upper;
00519             if (s<=0) { w += s; s = 0; b = c+lower; }
00520             else b = s*w+n;
00521             mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
00522             mrc.data = mrc.store+mrc.skip;
00523 
00524             if (+(mrc.cw*LoadOnEntry)) {
00525                 REPORT
00526                     Real* ColCopy = mrc.data; Real* Mstore = store+b;
00527                 while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
00528             }
00529     }
00530 
00531     void BandMatrix::RestoreCol(MatrixRowCol& mrc) {
00532         REPORT
00533             int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
00534             Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
00535             Real* Cstore = mrc.data;
00536             int w = mrc.storage;
00537             while (w--) { *Mstore = *Cstore++; Mstore += n; }
00538     }
00539 
00540     // routines for symmetric band matrix
00541 
00542     void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc) {
00543         REPORT
00544             int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
00545             mrc.length = ncols;
00546             if (s<0) { w1 += s; o -= s; s = 0; }
00547             mrc.skip = s;
00548 
00549             if (+(mrc.cw*DirectPart))
00550             { REPORT  mrc.data = store+o; mrc.storage = w1; }
00551             else {
00552                 // don't allow StoreOnExit and !DirectPart
00553                 if (+(mrc.cw*StoreOnExit))
00554                     Throw(InternalException("SymmetricBandMatrix::GetRow(MatrixRowCol&)"));
00555                 int w = w1+lower; s += w-ncols; Real* RowCopy;
00556                 if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
00557                 if (!(mrc.cw*HaveStore)) {
00558                     REPORT
00559                         RowCopy = new Real [2*lower+1]; MatrixErrorNoSpace(RowCopy);
00560                     MONITOR_REAL_NEW("Make (SmBGetRow)",2*lower+1,RowCopy)
00561                         mrc.cw += HaveStore; mrc.data = RowCopy;
00562                 }
00563                 else { REPORT  RowCopy = mrc.data; }
00564 
00565                 if (+(mrc.cw*LoadOnEntry)) {
00566                     REPORT
00567                         Real* Mstore = store+o;
00568                     while (w1--) *RowCopy++ = *Mstore++;
00569                     Mstore--;
00570                     while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
00571                 }
00572             }
00573     }
00574 
00575     void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc) {
00576         // don't allow StoreHere
00577         if (+(mrc.cw*StoreHere))
00578             Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
00579 
00580         int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows;
00581         REPORT
00582             int s = c-lower; int o = c*w1;
00583         if (s<0) { w1 += s; o -= s; s = 0; }
00584         mrc.skip = s;
00585 
00586         if (+(mrc.cw*DirectPart))
00587         { REPORT  mrc.data = store+o; mrc.storage = w1; }
00588         else {
00589             // don't allow StoreOnExit and !DirectPart
00590             if (+(mrc.cw*StoreOnExit))
00591                 Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
00592             int w = w1+lower; s += w-ncols; Real* ColCopy;
00593             if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
00594 
00595             if ( +(mrc.cw*HaveStore) ) { REPORT ColCopy = mrc.data; }
00596             else {
00597                 REPORT ColCopy = new Real [2*lower+1]; MatrixErrorNoSpace(ColCopy);
00598                 MONITOR_REAL_NEW("Make (SmBGetCol)",2*lower+1,ColCopy)
00599                     mrc.cw += HaveStore; mrc.data = ColCopy;
00600             }
00601 
00602             if (+(mrc.cw*LoadOnEntry)) {
00603                 REPORT
00604                     Real* Mstore = store+o;
00605                 while (w1--) *ColCopy++ = *Mstore++;
00606                 Mstore--;
00607                 while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
00608             }
00609         }
00610     }
00611 
00612     void SymmetricBandMatrix::GetCol(MatrixColX& mrc) {
00613         int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows;
00614         if (+(mrc.cw*DirectPart)) {
00615             REPORT
00616                 int b = c*w1+lower;
00617             mrc.skip = c; c += w1-nrows; w1 -= c; mrc.storage = w1;
00618             Real* ColCopy = mrc.data = mrc.store+mrc.skip;
00619             if (+(mrc.cw*LoadOnEntry)) {
00620                 REPORT
00621                     Real* Mstore = store+b;
00622                 while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
00623             }
00624         }
00625         else {
00626             REPORT
00627                 // don't allow StoreOnExit and !DirectPart
00628                 if (+(mrc.cw*StoreOnExit))
00629                     Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixColX&)"));
00630             int s = c-lower; int o = c*w1;
00631             if (s<0) { w1 += s; o -= s; s = 0; }
00632             mrc.skip = s;
00633 
00634             int w = w1+lower; s += w-ncols;
00635             if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
00636 
00637             Real* ColCopy = mrc.data = mrc.store+mrc.skip;
00638 
00639             if (+(mrc.cw*LoadOnEntry)) {
00640                 REPORT
00641                     Real* Mstore = store+o;
00642                 while (w1--) *ColCopy++ = *Mstore++;
00643                 Mstore--;
00644                 while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
00645             }
00646 
00647         }
00648     }
00649 
00650     void SymmetricBandMatrix::RestoreCol(MatrixColX& mrc) {
00651         REPORT
00652             int c = mrc.rowcol;
00653         Real* Mstore = store + c*lower+c+lower;
00654         Real* Cstore = mrc.data; int w = mrc.storage;
00655         while (w--) { *Mstore = *Cstore++; Mstore += lower; }
00656     }
00657 
00658     // *************************** destructors *******************************
00659 
00660     MatrixRowCol::~MatrixRowCol() {
00661         if (+(cw*HaveStore)) {
00662             // don't know length
00663             MONITOR_REAL_DELETE("Free    (RowCol)",-1,data)
00664                 delete [] data;
00665         }
00666     }
00667 
00668     MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }
00669 
00670     MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
00671 
00672     MatrixColX::~MatrixColX() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
00673 
00674 #ifdef use_namespace
00675 }
00676 #endif
Generated on Fri Aug 23 01:37:06 2002 for VirtualU by
1.2.17