Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmat3.cpp

Go to the documentation of this file.
00001 //$$ newmat3.cpp        Matrix get and restore rows and columns
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
00004 
00005 //#define WANT_STREAM
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 //#define MONITOR(what,storage,store)
00024 //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
00025 
00026 #define MONITOR(what,store,storage) {}
00027 
00028 
00029 // Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
00030 //
00031 // LoadOnEntry:
00032 //    Load data into MatrixRow or Col dummy array under GetRow or GetCol
00033 // StoreOnExit:
00034 //    Restore data to original matrix under RestoreRow or RestoreCol
00035 // DirectPart:
00036 //    Load or restore only part directly stored; must be set with StoreOnExit
00037 //    Still have decide how to handle this with symmetric
00038 // StoreHere:
00039 //    used in columns only - store data at supplied storage address;
00040 //    used for GetCol, NextCol & RestoreCol. No need to fill out zeros
00041 // HaveStore:
00042 //    dummy array has been assigned (internal use only).
00043 
00044 // For symmetric matrices, treat columns as rows unless StoreHere is set;
00045 // then stick to columns as this will give better performance for doing
00046 // inverses
00047 
00048 // How components are used:
00049 
00050 // Use rows wherever possible in preference to columns
00051 
00052 // Columns without StoreHere are used in in-exact transpose, sum column,
00053 // multiply a column vector, and maybe in future access to column,
00054 // additional multiply functions, add transpose
00055 
00056 // Columns with StoreHere are used in exact transpose (not symmetric matrices
00057 // or vectors, load only)
00058 
00059 // Columns with MatrixColX (Store to full column) are used in inverse and solve
00060 
00061 // Functions required for each matrix class
00062 
00063 // GetRow(MatrixRowCol& mrc)
00064 // GetCol(MatrixRowCol& mrc)
00065 // GetCol(MatrixColX& mrc)
00066 // RestoreRow(MatrixRowCol& mrc)
00067 // RestoreCol(MatrixRowCol& mrc)
00068 // RestoreCol(MatrixColX& mrc)
00069 // NextRow(MatrixRowCol& mrc)
00070 // NextCol(MatrixRowCol& mrc)
00071 // NextCol(MatrixColX& mrc)
00072 
00073 // The Restore routines assume StoreOnExit has already been checked
00074 // Defaults for the Next routines are given below
00075 // Assume cannot have both !DirectPart && StoreHere for MatrixRowCol routines
00076 
00077 
00078 // Default NextRow and NextCol:
00079 // will work as a default but need to override NextRow for efficiency
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                                                // 423
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                                                // 423
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 // routines for matrix
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) )      // ColumnVector
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          //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_value; }
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       //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_value; }
00156       if (i) for (;;)
00157           { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_value; }
00158    }
00159 }
00160 
00161 void Matrix::RestoreCol(MatrixRowCol& mrc)
00162 {
00163    // always check StoreOnExit before calling RestoreCol
00164    REPORT                                   // 429
00165    if (+(mrc.cw*HaveStore))
00166    {
00167       REPORT                                // 426
00168       Real* Mstore = store+mrc.rowcol; int i=nrows_value;
00169       Real* Cstore = mrc.data;
00170       // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_value; }
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    // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_value; }
00181    if (i) for (;;)
00182       { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_value; }
00183 }
00184 
00185 void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
00186 
00187 void Matrix::NextCol(MatrixRowCol& mrc)
00188 {
00189    REPORT                                        // 632
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          //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_value; }
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          // while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_value; }
00220          if (i) for (;;)
00221             { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_value; }
00222       }
00223    }
00224    else { REPORT mrc.cw -= StoreOnExit; }
00225 }
00226 
00227 // routines for diagonal matrix
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))              // should not happen
00241       Throw(InternalException("DiagonalMatrix::GetCol(MatrixRowCol&)"));
00242    else  { REPORT mrc.data=store+mrc.skip; }
00243                                                       // not accessed
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                   // 800
00256 
00257 void DiagonalMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
00258                         // not accessed
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 // routines for upper triangular matrix
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                                              // not accessed
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
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   // while (i--) { *Mstore = *Cstore++; Mstore += --j; }
00326   if (i) for (;;)
00327      { *Mstore = *Cstore++; if (!(--i)) break; Mstore += --j; }
00328 }
00329 
00330 void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
00331                   // 722
00332 
00333 // routines for lower triangular matrix
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                                            // not accessed
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
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    //while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
00389    if (i) for (;;)
00390       { *Mstore = *Cstore++; if (!(--i)) break; Mstore += ++col; }
00391 }
00392 
00393 void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
00394                            //712
00395 
00396 // routines for symmetric matrix
00397 
00398 void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
00399 {
00400    REPORT                                                //571
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       // do not allow StoreOnExit and !DirectPart
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                                         // 544
00421          Real* Mstore = store+(row*(row+1))/2; int i = row;
00422          while (i--) *RowCopy++ = *Mstore++;
00423          i = ncols_value-row;
00424          // while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
00425          if (i) for (;;)
00426             { *RowCopy++ = *Mstore; if (!(--i)) break; Mstore += ++row; }
00427       }
00428    }
00429 }
00430 
00431 void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
00432 {
00433    // do not allow StoreHere
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))    // actually get row ??
00441       { REPORT mrc.storage=col+1; mrc.data=store+(col*(col+1))/2; }
00442    else
00443    {
00444       // do not allow StoreOnExit and !DirectPart
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                                      // not accessed
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          // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
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                           // not accessed
00481          Real* ColCopy = mrc.data;
00482          Real* Mstore = store+(col*(col+3))/2;
00483          // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00484          if (i) for (;;)
00485             { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00486       }
00487    }
00488    else
00489    {
00490       REPORT
00491       // do not allow StoreOnExit and !DirectPart
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          // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00504          if (i) for (;;)
00505             { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00506       }
00507    }
00508 }
00509 
00510 // Do not need RestoreRow because we do not allow !DirectPart && StoreOnExit
00511 
00512 void SymmetricMatrix::RestoreCol(MatrixColX& mrc)
00513 {
00514    REPORT
00515    // Really do restore column
00516    int col=mrc.rowcol; Real* Cstore = mrc.data;
00517    Real* Mstore = store+(col*(col+3))/2; int i = nrows_value-col;
00518    // while (i--) { *Mstore = *Cstore++; Mstore+= ++col; }
00519    if (i) for (;;)
00520       { *Mstore = *Cstore++; if (!(--i)) break; Mstore+= ++col; }
00521 
00522 }
00523 
00524 // routines for row vector
00525 
00526 void RowVector::GetCol(MatrixRowCol& mrc)
00527 {
00528    REPORT
00529    // do not allow StoreHere
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); }           // not accessed
00564 
00565 
00566 // routines for band matrices
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 // should make special versions of this for upper and lower band matrices
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       // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
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       // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
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    // while (w--) { *Mstore = *Cstore++; Mstore += n; }
00644    if (w) for (;;)
00645       { *Mstore = *Cstore++; if (!(--w)) break; Mstore += n; }
00646 }
00647 
00648 // routines for symmetric band matrix
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       // do not allow StoreOnExit and !DirectPart
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    // do not allow StoreHere
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       // do not allow StoreOnExit and !DirectPart
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          // while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
00742          if (w1) for (;;)
00743             { *ColCopy++ = *Mstore; if (!(--w1)) break; Mstore += lower; }
00744       }
00745    }
00746    else
00747    {
00748       REPORT
00749       // do not allow StoreOnExit and !DirectPart
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    // while (w--) { *Mstore = *Cstore++; Mstore += lower; }
00780    if (w) for (;;)
00781       { *Mstore = *Cstore++; if (!(--w)) break; Mstore += lower; }
00782 }
00783 
00784 // routines for identity matrix
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))              // should not happen
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();            // must increase mrc.data so need IncrDiag
00817    int t1 = +(mrc.cw*LoadOnEntry);
00818    if (t1 && mrc.rowcol < ncols_value) { REPORT *(mrc.data)=*store; }
00819 }
00820 
00821 
00822 
00823 
00824 // *************************** destructors *******************************
00825 
00826 MatrixRowCol::~MatrixRowCol()
00827 {
00828    if (+(cw*HaveStore))
00829    {
00830       MONITOR_REAL_DELETE("Free    (RowCol)",-1,data)  // do not know length
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 

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3