newmat3.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 // Copyright (C) 1991,2,3,4: R B Davies
00008 
00009 //#define WANT_STREAM
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 //#define MONITOR(what,storage,store)
00028 //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
00029 
00030 #define MONITOR(what,store,storage) {}
00031 
00032 
00033 // Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
00034 //
00035 // LoadOnEntry:
00036 //    Load data into MatrixRow or Col dummy array under GetRow or GetCol
00037 // StoreOnExit:
00038 //    Restore data to original matrix under RestoreRow or RestoreCol
00039 // DirectPart:
00040 //    Load or restore only part directly stored; must be set with StoreOnExit
00041 //    Still have decide how to handle this with symmetric
00042 // StoreHere:
00043 //    used in columns only - store data at supplied storage address;
00044 //    used for GetCol, NextCol & RestoreCol. No need to fill out zeros
00045 // HaveStore:
00046 //    dummy array has been assigned (internal use only).
00047 
00048 // For symmetric matrices, treat columns as rows unless StoreHere is set;
00049 // then stick to columns as this will give better performance for doing
00050 // inverses
00051 
00052 // How components are used:
00053 
00054 // Use rows wherever possible in preference to columns
00055 
00056 // Columns without StoreHere are used in in-exact transpose, sum column,
00057 // multiply a column vector, and maybe in future access to column,
00058 // additional multiply functions, add transpose
00059 
00060 // Columns with StoreHere are used in exact transpose (not symmetric matrices
00061 // or vectors, load only)
00062 
00063 // Columns with MatrixColX (Store to full column) are used in inverse and solve
00064 
00065 // Functions required for each matrix class
00066 
00067 // GetRow(MatrixRowCol& mrc)
00068 // GetCol(MatrixRowCol& mrc)
00069 // GetCol(MatrixColX& mrc)
00070 // RestoreRow(MatrixRowCol& mrc)
00071 // RestoreCol(MatrixRowCol& mrc)
00072 // RestoreCol(MatrixColX& mrc)
00073 // NextRow(MatrixRowCol& mrc)
00074 // NextCol(MatrixRowCol& mrc)
00075 // NextCol(MatrixColX& mrc)
00076 
00077 // The Restore routines assume StoreOnExit has already been checked
00078 // Defaults for the Next routines are given below
00079 // Assume cannot have both !DirectPart && StoreHere for MatrixRowCol routines
00080 
00081 
00082 // Default NextRow and NextCol:
00083 // will work as a default but need to override NextRow for efficiency
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                                                // 423
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                                                // 423
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 // routines for matrix
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) )      // ColumnVector
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          //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
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       //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
00160       if (i) for (;;)
00161           { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
00162    }
00163 }
00164 
00165 void Matrix::RestoreCol(MatrixRowCol& mrc)
00166 {
00167    // always check StoreOnExit before calling RestoreCol
00168    REPORT                                   // 429
00169    if (+(mrc.cw*HaveStore))
00170    {
00171       REPORT                                // 426
00172       Real* Mstore = store+mrc.rowcol; int i=nrows_val;
00173       Real* Cstore = mrc.data;
00174       // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_val; }
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    // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_val; }
00185    if (i) for (;;)
00186       { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_val; }
00187 }
00188 
00189 void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
00190 
00191 void Matrix::NextCol(MatrixRowCol& mrc)
00192 {
00193    REPORT                                        // 632
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          //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
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          // while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
00224          if (i) for (;;)
00225             { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
00226       }
00227    }
00228    else { REPORT mrc.cw -= StoreOnExit; }
00229 }
00230 
00231 // routines for diagonal matrix
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))              // should not happen
00245       Throw(InternalException("DiagonalMatrix::GetCol(MatrixRowCol&)"));
00246    else  { REPORT mrc.data=store+mrc.skip; }
00247                                                       // not accessed
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                         // 800
00260 
00261 void DiagonalMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
00262                         // not accessed
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 // routines for upper triangular matrix
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                                              // not accessed
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
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   // while (i--) { *Mstore = *Cstore++; Mstore += --j; }
00330   if (i) for (;;)
00331      { *Mstore = *Cstore++; if (!(--i)) break; Mstore += --j; }
00332 }
00333 
00334 void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
00335                                                       // 722
00336 
00337 // routines for lower triangular matrix
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                                            // not accessed
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
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       // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
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    //while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
00393    if (i) for (;;)
00394       { *Mstore = *Cstore++; if (!(--i)) break; Mstore += ++col; }
00395 }
00396 
00397 void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
00398                                                          //712
00399 
00400 // routines for symmetric matrix
00401 
00402 void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
00403 {
00404    REPORT                                                //571
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       // do not allow StoreOnExit and !DirectPart
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                                         // 544
00425          Real* Mstore = store+(row*(row+1))/2; int i = row;
00426          while (i--) *RowCopy++ = *Mstore++;
00427          i = ncols_val-row;
00428          // while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
00429          if (i) for (;;)
00430             { *RowCopy++ = *Mstore; if (!(--i)) break; Mstore += ++row; }
00431       }
00432    }
00433 }
00434 
00435 void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
00436 {
00437    // do not allow StoreHere
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))    // actually get row ??
00445       { REPORT mrc.storage=col+1; mrc.data=store+(col*(col+1))/2; }
00446    else
00447    {
00448       // do not allow StoreOnExit and !DirectPart
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                                      // not accessed
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          // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
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                           // not accessed
00485          Real* ColCopy = mrc.data;
00486          Real* Mstore = store+(col*(col+3))/2;
00487          // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00488          if (i) for (;;)
00489             { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00490       }
00491    }
00492    else
00493    {
00494       REPORT
00495       // do not allow StoreOnExit and !DirectPart
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          // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
00508          if (i) for (;;)
00509             { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
00510       }
00511    }
00512 }
00513 
00514 // Do not need RestoreRow because we do not allow !DirectPart && StoreOnExit
00515 
00516 void SymmetricMatrix::RestoreCol(MatrixColX& mrc)
00517 {
00518    REPORT
00519    // Really do restore column
00520    int col=mrc.rowcol; Real* Cstore = mrc.data;
00521    Real* Mstore = store+(col*(col+3))/2; int i = nrows_val-col;
00522    // while (i--) { *Mstore = *Cstore++; Mstore+= ++col; }
00523    if (i) for (;;)
00524       { *Mstore = *Cstore++; if (!(--i)) break; Mstore+= ++col; }
00525 
00526 }
00527 
00528 // routines for row vector
00529 
00530 void RowVector::GetCol(MatrixRowCol& mrc)
00531 {
00532    REPORT
00533    // do not allow StoreHere
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); }           // not accessed
00568 
00569 
00570 // routines for band matrices
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 // should make special versions of this for upper and lower band matrices
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       // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
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       // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
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    // while (w--) { *Mstore = *Cstore++; Mstore += n; }
00648    if (w) for (;;)
00649       { *Mstore = *Cstore++; if (!(--w)) break; Mstore += n; }
00650 }
00651 
00652 // routines for symmetric band matrix
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       // do not allow StoreOnExit and !DirectPart
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    // do not allow StoreHere
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       // do not allow StoreOnExit and !DirectPart
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          // while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
00746          if (w1) for (;;)
00747             { *ColCopy++ = *Mstore; if (!(--w1)) break; Mstore += lower_val; }
00748       }
00749    }
00750    else
00751    {
00752       REPORT
00753       // do not allow StoreOnExit and !DirectPart
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    // while (w--) { *Mstore = *Cstore++; Mstore += lower_val; }
00784    if (w) for (;;)
00785       { *Mstore = *Cstore++; if (!(--w)) break; Mstore += lower_val; }
00786 }
00787 
00788 // routines for identity matrix
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))              // should not happen
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();            // must increase mrc.data so need IncrDiag
00821    int t1 = +(mrc.cw*LoadOnEntry);
00822    if (t1 && mrc.rowcol < ncols_val) { REPORT *(mrc.data)=*store; }
00823 }
00824 
00825 
00826 
00827 
00828 // *************************** destructors *******************************
00829 
00830 MatrixRowCol::~MatrixRowCol()
00831 {
00832    if (+(cw*HaveStore))
00833    {
00834       MONITOR_REAL_DELETE("Free    (RowCol)",-1,data)  // do not know length
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 


kni
Author(s): Martin Günther
autogenerated on Thu Jun 6 2019 21:42:34