home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / sources / misc / 4253 / newmat6.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  13.0 KB  |  581 lines
  1. //$$ newmat6.cxx            Operators, element access, submatrices
  2.  
  3. // Copyright (C) 1991,2,3: R B Davies
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmat.h"
  8. #include "newmatrc.h"
  9.  
  10.  
  11. //#define REPORT { static ExeCounter ExeCount(__LINE__,6); ExeCount++; }
  12.  
  13. #define REPORT {}
  14.  
  15. /*************************** general utilities *************************/
  16.  
  17. static int tristore(int n)                      // els in triangular matrix
  18. { return (n*(n+1))/2; }
  19.  
  20.  
  21. /****************************** operators *******************************/
  22.  
  23. Real& Matrix::operator()(int m, int n)
  24. {
  25.    if (m<=0 || m>nrows || n<=0 || n>ncols)
  26.       Throw(IndexException(m,n,*this));
  27.    return store[(m-1)*ncols+n-1];
  28. }
  29.  
  30. Real& SymmetricMatrix::operator()(int m, int n)
  31. {
  32.    if (m<=0 || n<=0 || m>nrows || n>ncols)
  33.       Throw(IndexException(m,n,*this));
  34.    if (m>=n) return store[tristore(m-1)+n-1];
  35.    else return store[tristore(n-1)+m-1];
  36. }
  37.  
  38. Real& UpperTriangularMatrix::operator()(int m, int n)
  39. {
  40.    if (m<=0 || n<m || n>ncols)
  41.       Throw(IndexException(m,n,*this));
  42.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  43. }
  44.  
  45. Real& LowerTriangularMatrix::operator()(int m, int n)
  46. {
  47.    if (n<=0 || m<n || m>nrows)
  48.       Throw(IndexException(m,n,*this));
  49.    return store[tristore(m-1)+n-1];
  50. }
  51.  
  52. Real& DiagonalMatrix::operator()(int m, int n)
  53. {
  54.    if (n<=0 || m!=n || m>nrows || n>ncols)
  55.       Throw(IndexException(m,n,*this));
  56.    return store[n-1];
  57. }
  58.  
  59. Real& DiagonalMatrix::operator()(int m)
  60. {
  61.    if (m<=0 || m>nrows) Throw(IndexException(m,*this));
  62.    return store[m-1];
  63. }
  64.  
  65. Real& ColumnVector::operator()(int m)
  66. {
  67.    if (m<=0 || m> nrows) Throw(IndexException(m,*this));
  68.    return store[m-1];
  69. }
  70.  
  71. Real& RowVector::operator()(int n)
  72. {
  73.    if (n<=0 || n> ncols) Throw(IndexException(n,*this));
  74.    return store[n-1];
  75. }
  76.  
  77. Real& BandMatrix::operator()(int m, int n)
  78. {
  79.    int w = upper+lower+1; int i = lower+n-m;
  80.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  81.       Throw(IndexException(m,n,*this));
  82.    return store[w*(m-1)+i];
  83. }
  84.  
  85. Real& UpperBandMatrix::operator()(int m, int n)
  86. {
  87.    int w = upper+1; int i = n-m;
  88.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  89.       Throw(IndexException(m,n,*this));
  90.    return store[w*(m-1)+i];
  91. }
  92.  
  93. Real& LowerBandMatrix::operator()(int m, int n)
  94. {
  95.    int w = lower+1; int i = lower+n-m;
  96.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  97.       Throw(IndexException(m,n,*this));
  98.    return store[w*(m-1)+i];
  99. }
  100.  
  101. Real& SymmetricBandMatrix::operator()(int m, int n)
  102. {
  103.    int w = lower+1;
  104.    if (m>=n)
  105.    {
  106.       int i = lower+n-m;
  107.       if ( m>nrows || n<=0 || i<0 )
  108.          Throw(IndexException(m,n,*this));
  109.       return store[w*(m-1)+i];
  110.    }
  111.    else
  112.    {
  113.       int i = lower+m-n;
  114.       if ( n>nrows || m<=0 || i<0 )
  115.          Throw(IndexException(m,n,*this));
  116.       return store[w*(n-1)+i];
  117.    }
  118. }
  119.  
  120. #ifndef __ZTC__
  121.  
  122. Real Matrix::operator()(int m, int n) const
  123. {
  124.    if (m<=0 || m>nrows || n<=0 || n>ncols)
  125.       Throw(IndexException(m,n,*this));
  126.    return store[(m-1)*ncols+n-1];
  127. }
  128.  
  129. Real SymmetricMatrix::operator()(int m, int n) const
  130. {
  131.    if (m<=0 || n<=0 || m>nrows || n>ncols)
  132.       Throw(IndexException(m,n,*this));
  133.    if (m>=n) return store[tristore(m-1)+n-1];
  134.    else return store[tristore(n-1)+m-1];
  135. }
  136.  
  137. Real UpperTriangularMatrix::operator()(int m, int n) const
  138. {
  139.    if (m<=0 || n<m || n>ncols)
  140.       Throw(IndexException(m,n,*this));
  141.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  142. }
  143.  
  144. Real LowerTriangularMatrix::operator()(int m, int n) const
  145. {
  146.    if (n<=0 || m<n || m>nrows)
  147.       Throw(IndexException(m,n,*this));
  148.    return store[tristore(m-1)+n-1];
  149. }
  150.  
  151. Real DiagonalMatrix::operator()(int m, int n) const
  152. {
  153.    if (n<=0 || m!=n || m>nrows || n>ncols)
  154.       Throw(IndexException(m,n,*this));
  155.    return store[n-1];
  156. }
  157.  
  158. Real DiagonalMatrix::operator()(int m) const
  159. {
  160.    if (m<=0 || m>nrows) Throw(IndexException(m,*this));
  161.    return store[m-1];
  162. }
  163.  
  164. Real ColumnVector::operator()(int m) const
  165. {
  166.    if (m<=0 || m> nrows) Throw(IndexException(m,*this));
  167.    return store[m-1];
  168. }
  169.  
  170. Real RowVector::operator()(int n) const
  171. {
  172.    if (n<=0 || n> ncols) Throw(IndexException(n,*this));
  173.    return store[n-1];
  174. }
  175.  
  176. Real BandMatrix::operator()(int m, int n) const
  177. {
  178.    int w = upper+lower+1; int i = lower+n-m;
  179.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  180.       Throw(IndexException(m,n,*this));
  181.    return store[w*(m-1)+i];
  182. }
  183.  
  184. Real SymmetricBandMatrix::operator()(int m, int n) const
  185. {
  186.    int w = lower+1;
  187.    if (m>=n)
  188.    {
  189.       int i = lower+n-m;
  190.       if ( m>nrows || n<=0 || i<0 )
  191.          Throw(IndexException(m,n,*this));
  192.       return store[w*(m-1)+i];
  193.    }
  194.    else
  195.    {
  196.       int i = lower+m-n;
  197.       if ( n>nrows || m<=0 || i<0 )
  198.          Throw(IndexException(m,n,*this));
  199.       return store[w*(n-1)+i];
  200.    }
  201. }
  202.  
  203. #endif
  204.  
  205. Real BaseMatrix::AsScalar() const
  206. {
  207.    REPORT
  208.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  209.    if (gm->nrows!=1 || gm->ncols!=1)
  210.    {
  211.       Try
  212.       {
  213.          Tracer tr("AsScalar");
  214.      Throw(ProgramException("Cannot convert to scalar", *gm));
  215.       }
  216.       CatchAll { gm->tDelete(); Throw(); }
  217.    }
  218.    Real x = *(gm->store); gm->tDelete(); return x;
  219. }
  220.  
  221. #ifdef TEMPS_DESTROYED_QUICKLY
  222.  
  223. AddedMatrix& BaseMatrix::operator+(const BaseMatrix& bm) const
  224. {
  225.    REPORT
  226.    AddedMatrix* x = new AddedMatrix(this, &bm);
  227.    MatrixErrorNoSpace(x);
  228.    return *x;
  229. }
  230.  
  231. MultipliedMatrix& BaseMatrix::operator*(const BaseMatrix& bm) const
  232. {
  233.    REPORT
  234.    MultipliedMatrix* x = new MultipliedMatrix(this, &bm);
  235.    MatrixErrorNoSpace(x);
  236.    return *x;
  237. }
  238.  
  239. //SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx) const
  240. SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx)
  241. {
  242.    REPORT
  243.    SolvedMatrix* x;
  244.    Try { x = new SolvedMatrix(bm, &bmx); MatrixErrorNoSpace(x); }
  245.    CatchAll { delete this; Throw(); }
  246.    delete this;                // since we are using bm rather than this
  247.    return *x;
  248. }
  249.  
  250. SubtractedMatrix& BaseMatrix::operator-(const BaseMatrix& bm) const
  251. {
  252.    REPORT
  253.    SubtractedMatrix* x = new SubtractedMatrix(this, &bm);
  254.    MatrixErrorNoSpace(x);
  255.    return *x;
  256. }
  257.  
  258. ShiftedMatrix& BaseMatrix::operator+(Real f) const
  259. {
  260.    REPORT
  261.    ShiftedMatrix* x = new ShiftedMatrix(this, f);
  262.    MatrixErrorNoSpace(x);
  263.    return *x;
  264. }
  265.  
  266. ScaledMatrix& BaseMatrix::operator*(Real f) const
  267. {
  268.    REPORT
  269.    ScaledMatrix* x = new ScaledMatrix(this, f);
  270.    MatrixErrorNoSpace(x);
  271.    return *x;
  272. }
  273.  
  274. ScaledMatrix& BaseMatrix::operator/(Real f) const
  275. {
  276.    REPORT
  277.    ScaledMatrix* x = new ScaledMatrix(this, 1.0/f);
  278.    MatrixErrorNoSpace(x);
  279.    return *x;
  280. }
  281.  
  282. ShiftedMatrix& BaseMatrix::operator-(Real f) const
  283. {
  284.    REPORT
  285.    ShiftedMatrix* x = new ShiftedMatrix(this, -f);
  286.    MatrixErrorNoSpace(x);
  287.    return *x;
  288. }
  289.  
  290. TransposedMatrix& BaseMatrix::t() const
  291. {
  292.    REPORT
  293.    TransposedMatrix* x = new TransposedMatrix(this);
  294.    MatrixErrorNoSpace(x);
  295.    return *x;
  296. }
  297.  
  298. NegatedMatrix& BaseMatrix::operator-() const
  299. {
  300.    REPORT
  301.    NegatedMatrix* x = new NegatedMatrix(this);
  302.    MatrixErrorNoSpace(x);
  303.    return *x;
  304. }
  305.  
  306. InvertedMatrix& BaseMatrix::i() const
  307. {
  308.    REPORT
  309.    InvertedMatrix* x = new InvertedMatrix(this);
  310.    MatrixErrorNoSpace(x);
  311.    return *x;
  312. }
  313.  
  314. ConstMatrix& GeneralMatrix::c() const
  315. {
  316.    if (tag != -1)
  317.       Throw(ProgramException(".c() applied to temporary matrix"));
  318.    REPORT
  319.    ConstMatrix* x = new ConstMatrix(this);
  320.    MatrixErrorNoSpace(x);
  321.    return *x;
  322. }
  323.  
  324. RowedMatrix& BaseMatrix::AsRow() const
  325. {
  326.    REPORT
  327.    RowedMatrix* x = new RowedMatrix(this);
  328.    MatrixErrorNoSpace(x);
  329.    return *x;
  330. }
  331.  
  332. ColedMatrix& BaseMatrix::AsColumn() const
  333. {
  334.    REPORT
  335.    ColedMatrix* x = new ColedMatrix(this);
  336.    MatrixErrorNoSpace(x);
  337.    return *x;
  338. }
  339.  
  340. DiagedMatrix& BaseMatrix::AsDiagonal() const
  341. {
  342.    REPORT
  343.    DiagedMatrix* x = new DiagedMatrix(this);
  344.    MatrixErrorNoSpace(x);
  345.    return *x;
  346. }
  347.  
  348. MatedMatrix& BaseMatrix::AsMatrix(int nrx, int ncx) const
  349. {
  350.    REPORT
  351.    MatedMatrix* x = new MatedMatrix(this,nrx,ncx);
  352.    MatrixErrorNoSpace(x);
  353.    return *x;
  354. }
  355.  
  356. #else
  357.  
  358. AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
  359. { REPORT return AddedMatrix(this, &bm); }
  360.  
  361. MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
  362. { REPORT return MultipliedMatrix(this, &bm); }
  363.  
  364. SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
  365. { REPORT return SolvedMatrix(bm, &bmx); }
  366.  
  367. SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
  368. { REPORT return SubtractedMatrix(this, &bm); }
  369.  
  370. ShiftedMatrix BaseMatrix::operator+(Real f) const
  371. { REPORT return ShiftedMatrix(this, f); }
  372.  
  373. ScaledMatrix BaseMatrix::operator*(Real f) const
  374. { REPORT return ScaledMatrix(this, f); }
  375.  
  376. ScaledMatrix BaseMatrix::operator/(Real f) const
  377. { REPORT return ScaledMatrix(this, 1.0/f); }
  378.  
  379. ShiftedMatrix BaseMatrix::operator-(Real f) const
  380. { REPORT return ShiftedMatrix(this, -f); }
  381.  
  382. TransposedMatrix BaseMatrix::t() const
  383. { REPORT return TransposedMatrix(this); }
  384.  
  385. NegatedMatrix BaseMatrix::operator-() const
  386. { REPORT return NegatedMatrix(this); }
  387.  
  388. InvertedMatrix BaseMatrix::i() const
  389. { REPORT return InvertedMatrix(this); }
  390.  
  391. ConstMatrix GeneralMatrix::c() const
  392. {
  393.    if (tag != -1)
  394.       Throw(ProgramException(".c() applied to temporary matrix"));
  395.    REPORT return ConstMatrix(this);
  396. }
  397.  
  398. RowedMatrix BaseMatrix::AsRow() const
  399. { REPORT return RowedMatrix(this); }
  400.  
  401. ColedMatrix BaseMatrix::AsColumn() const
  402. { REPORT return ColedMatrix(this); }
  403.  
  404. DiagedMatrix BaseMatrix::AsDiagonal() const
  405. { REPORT return DiagedMatrix(this); }
  406.  
  407. MatedMatrix BaseMatrix::AsMatrix(int nrx, int ncx) const
  408. { REPORT return MatedMatrix(this,nrx,ncx); }
  409.  
  410. #endif
  411.  
  412. void GeneralMatrix::operator=(Real f)
  413. { REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }
  414.  
  415. void Matrix::operator=(const BaseMatrix& X)
  416. {
  417.    REPORT //CheckConversion(X);
  418.    MatrixConversionCheck mcc;
  419.    Eq(X,MatrixType::Rt);
  420. } 
  421.  
  422. void RowVector::operator=(const BaseMatrix& X)
  423. {
  424.    REPORT  // CheckConversion(X);
  425.    MatrixConversionCheck mcc;
  426.    Eq(X,MatrixType::RV);
  427.    if (nrows!=1)
  428.    {
  429.       Tracer tr("RowVector(=)");
  430.       Throw(VectorException(*this));
  431.    }
  432. }
  433.  
  434. void ColumnVector::operator=(const BaseMatrix& X)
  435. {
  436.    REPORT //CheckConversion(X);
  437.    MatrixConversionCheck mcc;
  438.    Eq(X,MatrixType::CV);
  439.    if (ncols!=1)
  440.    {
  441.       Tracer tr("ColumnVector(=)");
  442.       Throw(VectorException(*this));
  443.    }
  444. }
  445.  
  446. void SymmetricMatrix::operator=(const BaseMatrix& X)
  447. {
  448.    REPORT // CheckConversion(X);
  449.    MatrixConversionCheck mcc;
  450.    Eq(X,MatrixType::Sm);
  451. }
  452.  
  453. void UpperTriangularMatrix::operator=(const BaseMatrix& X)
  454. {
  455.    REPORT //CheckConversion(X);
  456.    MatrixConversionCheck mcc;
  457.    Eq(X,MatrixType::UT);
  458. }
  459.  
  460. void LowerTriangularMatrix::operator=(const BaseMatrix& X)
  461. {
  462.    REPORT //CheckConversion(X);
  463.    MatrixConversionCheck mcc;
  464.    Eq(X,MatrixType::LT);
  465. }
  466.  
  467. void DiagonalMatrix::operator=(const BaseMatrix& X)
  468. {
  469.    REPORT // CheckConversion(X);
  470.    MatrixConversionCheck mcc;
  471.    Eq(X,MatrixType::Dg);
  472. }
  473.  
  474. void GeneralMatrix::operator<<(const Real* r)
  475. {
  476.    REPORT
  477.    int i = storage; Real* s=store;
  478.    while(i--) *s++ = *r++;
  479. }
  480.  
  481.  
  482. /************************* element access *********************************/
  483.  
  484. Real& Matrix::element(int m, int n)
  485. {
  486.    if (m<0 || m>= nrows || n<0 || n>= ncols)
  487.       Throw(IndexException(m,n,*this,TRUE));
  488.    return store[m*ncols+n];
  489. }
  490.  
  491. Real& SymmetricMatrix::element(int m, int n)
  492. {
  493.    if (m<0 || n<0 || m >= nrows || n>=ncols)
  494.       Throw(IndexException(m,n,*this,TRUE));
  495.    if (m>=n) return store[tristore(m)+n];
  496.    else return store[tristore(n)+m];
  497. }
  498.  
  499. Real& UpperTriangularMatrix::element(int m, int n)
  500. {
  501.    if (m<0 || n<m || n>=ncols)
  502.       Throw(IndexException(m,n,*this,TRUE));
  503.    return store[m*ncols+n-tristore(m)];
  504. }
  505.  
  506. Real& LowerTriangularMatrix::element(int m, int n)
  507. {
  508.    if (n<0 || m<n || m>=nrows)
  509.       Throw(IndexException(m,n,*this,TRUE));
  510.    return store[tristore(m)+n];
  511. }
  512.  
  513. Real& DiagonalMatrix::element(int m, int n)
  514. {
  515.    if (n<0 || m!=n || m>=nrows || n>=ncols)
  516.       Throw(IndexException(m,n,*this,TRUE));
  517.    return store[n];
  518. }
  519.  
  520. Real& DiagonalMatrix::element(int m)
  521. {
  522.    if (m<0 || m>=nrows) Throw(IndexException(m,*this,TRUE));
  523.    return store[m];
  524. }
  525.  
  526. Real& ColumnVector::element(int m)
  527. {
  528.    if (m<0 || m>= nrows) Throw(IndexException(m,*this,TRUE));
  529.    return store[m];
  530. }
  531.  
  532. Real& RowVector::element(int n)
  533. {
  534.    if (n<0 || n>= ncols)  Throw(IndexException(n,*this,TRUE));
  535.    return store[n];
  536. }
  537.  
  538. Real& BandMatrix::element(int m, int n)
  539. {
  540.    int w = upper+lower+1; int i = lower+n-m;
  541.    if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  542.       Throw(IndexException(m,n,*this,TRUE));
  543.    return store[w*m+i];
  544. }
  545.  
  546. Real& UpperBandMatrix::element(int m, int n)
  547. {
  548.    int w = upper+1; int i = n-m;
  549.    if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  550.       Throw(IndexException(m,n,*this,TRUE));
  551.    return store[w*m+i];
  552. }
  553.  
  554. Real& LowerBandMatrix::element(int m, int n)
  555. {
  556.    int w = lower+1; int i = lower+n-m;
  557.    if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  558.       Throw(IndexException(m,n,*this,TRUE));
  559.    return store[w*m+i];
  560. }
  561.  
  562. Real& SymmetricBandMatrix::element(int m, int n)
  563. {
  564.    int w = lower+1;
  565.    if (m>=n)
  566.    {
  567.       int i = lower+n-m;
  568.       if ( m>=nrows || n<0 || i<0 )
  569.          Throw(IndexException(m,n,*this,TRUE));
  570.       return store[w*m+i];
  571.    }
  572.    else
  573.    {
  574.       int i = lower+m-n;
  575.       if ( n>=nrows || m<0 || i<0 )
  576.          Throw(IndexException(m,n,*this,TRUE));
  577.       return store[w*n+i];
  578.    }
  579. }
  580.  
  581.