home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / NEWMAT5.CXX < prev    next >
C/C++ Source or Header  |  1992-11-29  |  11KB  |  386 lines
  1. //$$ newmat5.cxx         Transpose, evaluate etc
  2.  
  3. // Copyright (C) 1991,2: R B Davies
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmat.h"
  8. #include "newmatrc.h"
  9.  
  10. //#define REPORT { static ExeCounter ExeCount(__LINE__,5); ExeCount++; }
  11.  
  12. #define REPORT {}
  13.  
  14.  
  15. /************************ carry out operations ******************************/
  16.  
  17.  
  18. GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
  19. {
  20.    if (!mt) mt = Type().t();           // type of transposed matrix
  21.    GeneralMatrix* gm1 = mt.New(ncols,nrows,tm); gm1->ReleaseAndDelete();
  22.  
  23.    if (mt == Type().t())
  24.    {
  25.       REPORT
  26.       for (int i=0; i<ncols; i++)
  27.       {
  28.      MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
  29.          MatrixCol mc(this, mr.Store(), LoadOnEntry, i);
  30.       }
  31.    }
  32.    else
  33.    {
  34.       REPORT
  35.       MatrixRow mr(this, LoadOnEntry);
  36.       MatrixCol mc(gm1, StoreOnExit+DirectPart);
  37.       int i = nrows;
  38.       while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
  39.    }
  40.    tDelete(); return gm1;
  41. }
  42.  
  43. GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  44. { REPORT  return Evaluate(mt); }
  45.  
  46.  
  47. GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  48. { REPORT return Evaluate(mt); }
  49.  
  50. Boolean GeneralMatrix::IsZero() const
  51. {
  52.    REPORT
  53.    Real* s=store; int i=storage;
  54.    while (i--) { if (*s++) return FALSE; }
  55.    return TRUE;
  56. }
  57.  
  58. GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
  59. {
  60.    REPORT
  61.    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  62.    gmx->nrows = 1; gmx->ncols = gmx->storage = storage;
  63.    return BorrowStore(gmx,mt);
  64. }
  65.  
  66. GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
  67. {
  68.    REPORT
  69.    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  70.    gmx->ncols = 1; gmx->nrows = gmx->storage = storage;
  71.    return BorrowStore(gmx,mt);
  72. }
  73.  
  74. GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
  75. {
  76.    if (!mt) { REPORT return this; }
  77.    if (mt==this->Type()) { REPORT return this; }
  78.    REPORT
  79.    GeneralMatrix* gmx = mt.New(nrows,ncols,this);
  80.    MatrixRow mr(this, LoadOnEntry);
  81.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  82.    int i=nrows;
  83.    while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  84.    tDelete(); gmx->ReleaseAndDelete(); return gmx;
  85. }
  86.  
  87. GeneralMatrix* ConstMatrix::Evaluate(MatrixType mt)
  88. {
  89.    if (!mt || mt==cgm->Type())
  90.    {
  91.       REPORT
  92. #ifdef TEMPS_DESTROYED_QUICKLY
  93.       GeneralMatrix* gmx = (GeneralMatrix*)cgm; delete this; return gmx;
  94. #else
  95.       return (GeneralMatrix*)cgm;
  96. #endif
  97.    }
  98.    REPORT
  99.    GeneralMatrix* gmx = cgm->Type().New(cgm->Nrows(),cgm->Ncols(),this);
  100.    MatrixRow mr((GeneralMatrix*)cgm, LoadOnEntry);//assume won't change this
  101.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  102.    int i=cgm->Nrows();
  103.    while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  104.    gmx->ReleaseAndDelete();
  105. #ifdef TEMPS_DESTROYED_QUICKLY
  106.    delete this;
  107. #endif
  108.    return gmx; // no tDelete
  109. }
  110.  
  111. GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
  112. {
  113.    gm=((BaseMatrix*&)bm)->Evaluate();
  114.    if (!mt) mt = gm->Type().AddEqualEl();
  115.    int nr=gm->Nrows(); int nc=gm->Ncols();
  116.    if (!(mt==gm->Type()))
  117.    {
  118.       REPORT
  119.       GeneralMatrix* gmx = mt.New(nr,nc,this);
  120.       MatrixRow mr(gm, LoadOnEntry); 
  121.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  122.       while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
  123.       gmx->ReleaseAndDelete(); gm->tDelete();
  124. #ifdef TEMPS_DESTROYED_QUICKLY
  125.       delete this;
  126. #endif
  127.       return gmx;
  128.    }
  129.    else if (gm->reuse())
  130.    {
  131.       REPORT gm->Add(f);
  132. #ifdef TEMPS_DESTROYED_QUICKLY
  133.       GeneralMatrix* gmx = gm; delete this; return gmx;
  134. #else
  135.       return gm;
  136. #endif
  137.    }
  138.    else
  139.    {
  140.       REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
  141.       gmy->ReleaseAndDelete(); gmy->Add(gm,f);
  142. #ifdef TEMPS_DESTROYED_QUICKLY
  143.       delete this;
  144. #endif
  145.       return gmy;
  146.    }
  147. }
  148.  
  149. GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
  150. {
  151.    gm=((BaseMatrix*&)bm)->Evaluate();
  152.    int nr=gm->Nrows(); int nc=gm->Ncols();
  153.    if (!mt) mt = gm->Type();
  154.    if (mt==gm->Type())
  155.    {
  156.       if (gm->reuse())
  157.       {
  158.          REPORT gm->Multiply(f);
  159. #ifdef TEMPS_DESTROYED_QUICKLY
  160.          GeneralMatrix* gmx = gm; delete this; return gmx;
  161. #else
  162.          return gm;
  163. #endif
  164.       }
  165.       else
  166.       {
  167.          REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
  168.          gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
  169. #ifdef TEMPS_DESTROYED_QUICKLY
  170.          delete this;
  171. #endif
  172.          return gmx;
  173.       }
  174.    }
  175.    else
  176.    {
  177.       REPORT
  178.       GeneralMatrix* gmx = mt.New(nr,nc,this);
  179.       MatrixRow mr(gm, LoadOnEntry); 
  180.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  181.       while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
  182.       gmx->ReleaseAndDelete(); gm->tDelete();
  183. #ifdef TEMPS_DESTROYED_QUICKLY
  184.       delete this;
  185. #endif
  186.       return gmx;
  187.    }
  188. }
  189.  
  190. GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
  191. {
  192.    gm=((BaseMatrix*&)bm)->Evaluate();
  193.    if (!mt) mt = gm->Type();
  194.    int nr=gm->Nrows(); int nc=gm->Ncols();
  195.    if (mt==gm->Type())
  196.    {
  197.       if (gm->reuse())
  198.       {
  199.          REPORT gm->Negate();
  200. #ifdef TEMPS_DESTROYED_QUICKLY
  201.          GeneralMatrix* gmx = gm; delete this; return gmx;
  202. #else
  203.          return gm;
  204. #endif
  205.       }
  206.       else
  207.       {
  208.          REPORT
  209.          GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
  210.          gmx->ReleaseAndDelete(); gmx->Negate(gm);
  211. #ifdef TEMPS_DESTROYED_QUICKLY
  212.          delete this;
  213. #endif
  214.          return gmx;
  215.       }
  216.    }
  217.    else
  218.    {
  219.       REPORT
  220.       GeneralMatrix* gmx = mt.New(nr,nc,this);
  221.       MatrixRow mr(gm, LoadOnEntry); 
  222.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  223.       while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
  224.       gmx->ReleaseAndDelete(); gm->tDelete();
  225. #ifdef TEMPS_DESTROYED_QUICKLY
  226.       delete this;
  227. #endif
  228.       return gmx;
  229.    }
  230. }   
  231.  
  232. GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
  233. {
  234.    REPORT
  235.    gm=((BaseMatrix*&)bm)->Evaluate();
  236.    if (!mt) mt = gm->Type().t();
  237.    GeneralMatrix* gmx=gm->Transpose(this, mt);
  238. #ifdef TEMPS_DESTROYED_QUICKLY
  239.    delete this;
  240. #endif
  241.    return gmx;
  242. }
  243.    
  244. GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
  245. {
  246.    gm = ((BaseMatrix*&)bm)->Evaluate();
  247.    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  248.    gmx->nrows = 1; gmx->ncols = gmx->storage = gm->storage;
  249. #ifdef TEMPS_DESTROYED_QUICKLY
  250.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  251. #else
  252.    return gm->BorrowStore(gmx,mt);
  253. #endif
  254. }
  255.  
  256. GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
  257. {
  258.    gm = ((BaseMatrix*&)bm)->Evaluate();
  259.    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  260.    gmx->ncols = 1; gmx->nrows = gmx->storage = gm->storage;
  261. #ifdef TEMPS_DESTROYED_QUICKLY
  262.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  263. #else
  264.    return gm->BorrowStore(gmx,mt);
  265. #endif
  266. }
  267.  
  268. GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
  269. {
  270.    gm = ((BaseMatrix*&)bm)->Evaluate();
  271.    GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
  272.    gmx->nrows = gmx->ncols = gmx->storage = gm->storage;
  273. #ifdef TEMPS_DESTROYED_QUICKLY
  274.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  275. #else
  276.    return gm->BorrowStore(gmx,mt);
  277. #endif
  278. }
  279.  
  280. GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
  281. {
  282.    Tracer tr("MatedMatrix::Evaluate");
  283.    gm = ((BaseMatrix*&)bm)->Evaluate();
  284.    GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
  285.    gmx->nrows = nr; gmx->ncols = nc; gmx->storage = gm->storage;
  286.    if (nr*nc != gmx->storage)
  287.       Throw(IncompatibleDimensionsException());
  288. #ifdef TEMPS_DESTROYED_QUICKLY
  289.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  290. #else
  291.    return gm->BorrowStore(gmx,mt);
  292. #endif
  293. }
  294.  
  295. GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
  296. {
  297.    REPORT
  298.    Tracer tr("SubMatrix(evaluate)");
  299.    gm = ((BaseMatrix*&)bm)->Evaluate();
  300.    if (row_number < 0) row_number = gm->Nrows();
  301.    if (col_number < 0) col_number = gm->Ncols();
  302.    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  303.       Throw(SubMatrixDimensionException());
  304.    if (!mt) mt = MatrixType::Rt;
  305.    GeneralMatrix* gmx = mt.New(row_number, col_number, this);
  306.    int i = row_number;
  307.    MatrixRow mr(gm, LoadOnEntry, row_skip); 
  308.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  309.    MatrixRowCol sub;
  310.    while (i--)
  311.    {
  312.       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  313.       mrx.Copy(sub); mrx.Next(); mr.Next();
  314.    }
  315.    gmx->ReleaseAndDelete(); gm->tDelete();
  316. #ifdef TEMPS_DESTROYED_QUICKLY
  317.    delete this;
  318. #endif
  319.    return gmx;
  320. }   
  321.  
  322.  
  323. GeneralMatrix* ReturnMatrixX::Evaluate(MatrixType mt)
  324. {
  325. #ifdef TEMPS_DESTROYED_QUICKLY
  326.    GeneralMatrix* gmx = gm; delete this; return gmx->Evaluate(mt);
  327. #else
  328.    return gm->Evaluate(mt);
  329. #endif
  330. }
  331.  
  332.  
  333. void GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
  334. {
  335.    REPORT
  336.    Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
  337.    while (i--)
  338.    { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
  339.    i = storage & 3; while (i--) *s++ = *s1++ + f;
  340. }
  341.    
  342. void GeneralMatrix::Add(Real f)
  343. {
  344.    REPORT
  345.    Real* s=store; int i=(storage >> 2);
  346.    while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
  347.    i = storage & 3; while (i--) *s++ += f;
  348. }
  349.    
  350. void GeneralMatrix::Negate(GeneralMatrix* gm1)
  351. {
  352.    // change sign of elements
  353.    REPORT
  354.    Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
  355.    while (i--)
  356.    { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
  357.    i = storage & 3; while(i--) *s++ = -(*s1++);
  358. }
  359.    
  360. void GeneralMatrix::Negate()
  361. {
  362.    REPORT
  363.    Real* s=store; int i=(storage >> 2);
  364.    while (i--)
  365.    { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
  366.    i = storage & 3; while(i--) { *s = -(*s); s++; }
  367. }
  368.    
  369. void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
  370. {
  371.    REPORT
  372.    Real* s1=gm1->store; Real* s=store;  int i=(storage >> 2);
  373.    while (i--)
  374.    { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
  375.    i = storage & 3; while (i--) *s++ = *s1++ * f;
  376. }
  377.    
  378. void GeneralMatrix::Multiply(Real f)
  379. {
  380.    REPORT
  381.    Real* s=store; int i=(storage >> 2);
  382.    while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
  383.    i = storage & 3; while (i--) *s++ *= f;
  384. }
  385.    
  386.