home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / NEWMAT7.CXX < prev    next >
C/C++ Source or Header  |  1992-12-02  |  15KB  |  497 lines
  1. //$$ newmat7.cxx     Invert, solve, binary operations
  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__,7); ExeCount++; }
  11.  
  12. #define REPORT {}
  13.  
  14.  
  15. /***************************** solve routines ******************************/
  16.  
  17. GeneralMatrix* GeneralMatrix::MakeSolver()
  18. {
  19.    REPORT
  20.    GeneralMatrix* gm = new CroutMatrix(*this);
  21.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  22. }
  23.  
  24. GeneralMatrix* Matrix::MakeSolver()
  25. {
  26.    REPORT
  27.    GeneralMatrix* gm = new CroutMatrix(*this);
  28.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  29. }
  30.  
  31. void CroutMatrix::Solver(MatrixRowCol& mcout, const MatrixRowCol& mcin)
  32. {
  33.    REPORT
  34.    Real* el = mcin.store; int i = mcin.skip;
  35.    while (i--) *el++ = 0.0;
  36.    el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
  37.    while (i--) *el++ = 0.0;
  38.    lubksb(mcin.store, mcout.skip);
  39. }
  40.  
  41.  
  42. // Do we need check for entirely zero output?
  43.  
  44. void UpperTriangularMatrix::Solver(MatrixRowCol& mcout,
  45.    const MatrixRowCol& mcin)
  46. {
  47.    REPORT
  48.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  49.    while (i-- > 0) *elx++ = 0.0;
  50.    int nr = mcin.skip+mcin.storage; elx = mcin.store+nr; Real* el = elx;
  51.    int j = mcout.skip+mcout.storage-nr; int nc = ncols-nr; i = nr-mcout.skip;
  52.    while (j-- > 0) *elx++ = 0.0;
  53.    Real* Ael = store + (nr*(2*ncols-nr+1))/2; j = 0;
  54.    while (i-- > 0)
  55.    {
  56.       elx = el; Real sum = 0.0; int jx = j++; Ael -= nc;
  57.       while (jx--) sum += *(--Ael) * *(--elx);
  58.       elx--; *elx = (*elx - sum) / *(--Ael);
  59.    }
  60. }
  61.  
  62. void LowerTriangularMatrix::Solver(MatrixRowCol& mcout,
  63.    const MatrixRowCol& mcin)
  64. {
  65.    REPORT
  66.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  67.    while (i-- > 0) *elx++ = 0.0;
  68.    int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.store+i;
  69.    int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
  70.    while (j-- > 0) *elx++ = 0.0;
  71.    Real* el = mcin.store+nc; Real* Ael = store + (nc*(nc+1))/2; j = 0;
  72.    while (i-- > 0)
  73.    {
  74.       elx = el; Real sum = 0.0; int jx = j++; Ael += nc;
  75.       while (jx--) sum += *Ael++ * *elx++;
  76.       *elx = (*elx - sum) / *Ael++;
  77.    }
  78. }
  79.  
  80. /******************* carry out binary operations *************************/
  81.  
  82. static GeneralMatrix*
  83.    GeneralAdd(GeneralMatrix*,GeneralMatrix*,AddedMatrix*,MatrixType);
  84. static GeneralMatrix*
  85.    GeneralSub(GeneralMatrix*,GeneralMatrix*,SubtractedMatrix*,MatrixType);
  86. static GeneralMatrix*
  87.    GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
  88. static GeneralMatrix*
  89.    GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
  90.  
  91. GeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
  92. {
  93.    REPORT
  94.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  95.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  96. #ifdef TEMPS_DESTROYED_QUICKLY
  97.    GeneralMatrix* gmx = GeneralAdd(gm1,gm2,this,mt); delete this; return gmx;
  98. #else
  99.    return GeneralAdd(gm1,gm2,this,mt);
  100. #endif   
  101. }
  102.  
  103. GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
  104. {
  105.    REPORT
  106.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  107.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  108. #ifdef TEMPS_DESTROYED_QUICKLY
  109.    GeneralMatrix* gmx = GeneralSub(gm1,gm2,this,mt); delete this; return gmx;
  110. #else
  111.    return GeneralSub(gm1,gm2,this,mt);
  112. #endif   
  113. }
  114.  
  115. GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
  116. {
  117.    REPORT
  118. //   GeneralMatrix*
  119.    gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  120.    gm2 = gm2->Evaluate(gm2->Type().MultRHS());     // no symmetric on RHS
  121.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  122. #ifdef TEMPS_DESTROYED_QUICKLY
  123.    GeneralMatrix* gmx = GeneralMult(gm1, gm2, this, mt);
  124.    delete this; return gmx;
  125. #else
  126.    return GeneralMult(gm1, gm2, this, mt);
  127. #endif   
  128. }
  129.  
  130. GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
  131. {
  132.    REPORT
  133.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  134.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  135. #ifdef TEMPS_DESTROYED_QUICKLY
  136.    GeneralMatrix* gmx;
  137.    Try { gmx = GeneralSolv(gm1,gm2,this,mt); }
  138.    CatchAll { delete this; Throw(); }
  139.    delete this; return gmx;
  140. #else
  141.    return GeneralSolv(gm1,gm2,this,mt);
  142. #endif   
  143. }
  144.  
  145. // routines for adding or subtracting matrices of identical storage structure
  146.  
  147. static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  148. {
  149.    REPORT
  150.    Real* s1=gm1->Store(); Real* s2=gm2->Store();
  151.    Real* s=gm->Store(); int i=gm->Storage() >> 2;
  152.    while (i--)
  153.    {
  154.        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  155.        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  156.    }
  157.    i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
  158. }
  159.    
  160. static void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
  161. {
  162.    REPORT
  163.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  164.    while (i--)
  165.    { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
  166.    i=gm->Storage() & 3; while (i--) *s++ += *s2++;
  167. }
  168.  
  169. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  170. {
  171.    REPORT
  172.    Real* s1=gm1->Store(); Real* s2=gm2->Store();
  173.    Real* s=gm->Store(); int i=gm->Storage() >> 2;
  174.    while (i--)
  175.    {
  176.        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  177.        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  178.    }
  179.    i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
  180. }
  181.  
  182. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  183. {
  184.    REPORT
  185.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  186.    while (i--)
  187.    { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
  188.    i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
  189. }
  190.  
  191. static void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  192. {
  193.    REPORT
  194.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  195.    while (i--)
  196.    {
  197.       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  198.       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  199.    }
  200.    i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
  201. }
  202.  
  203. // routines for adding or subtracting matrices of different storage structure
  204.  
  205. static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  206. {
  207.    int nr = gm->Nrows();
  208.    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  209.    MatrixRow mr(gm, StoreOnExit+DirectPart);
  210.    while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  211. }
  212.    
  213. static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  214. // Add into first argument
  215. {
  216.    int nr = gm->Nrows();
  217.    MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
  218.    MatrixRow mr2(gm2, LoadOnEntry);
  219.    while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
  220. }
  221.  
  222. static void SubtractDS
  223.    (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  224. {
  225.    int nr = gm->Nrows();
  226.    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  227.    MatrixRow mr(gm, StoreOnExit+DirectPart);
  228.    while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  229. }
  230.  
  231. static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  232. {
  233.    int nr = gm->Nrows();
  234.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  235.    MatrixRow mr2(gm2, LoadOnEntry);
  236.    while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
  237. }
  238.  
  239. static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  240. {
  241.    int nr = gm->Nrows();
  242.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  243.    MatrixRow mr2(gm2, LoadOnEntry);
  244.    while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
  245. }
  246.  
  247. #ifdef GXX
  248. void AddedMatrix::SelectVersion
  249.    (MatrixType mtx, int& c1, int& c2) const
  250. #else
  251. void AddedMatrix::SelectVersion
  252.    (MatrixType mtx, Boolean& c1, Boolean& c2) const
  253. #endif
  254. // for determining version of add and subtract routines
  255. // will need to modify if further matrix structures are introduced
  256. {
  257.    MatrixBandWidth bm1 = gm1->BandWidth();
  258.    MatrixBandWidth bm2 = gm2->BandWidth();
  259.    MatrixBandWidth bmx = bm1 + bm2;
  260.    c1 = (mtx == gm1->Type()) && (bmx == bm1);
  261.    c2 = (mtx == gm2->Type()) && (bmx == bm2);
  262. }
  263.  
  264. static GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
  265.    AddedMatrix* am, MatrixType mtx)
  266. {
  267.    Tracer tr("GeneralAdd");
  268.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  269.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols()) 
  270.       Throw(IncompatibleDimensionsException());
  271.    if (!mtx) mtx = gm1->Type() + gm2->Type();
  272. #ifdef GXX
  273.    int c1,c2; am->SelectVersion(mtx,c1,c2);
  274. #else
  275.    Boolean c1,c2; am->SelectVersion(mtx,c1,c2); // causes problems for g++
  276. #endif
  277.    if (c1 && c2)
  278.    {
  279.       if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
  280.       else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
  281.       else
  282.       {
  283.          REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
  284.          gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
  285.       }
  286.    }
  287.    else
  288.    {
  289.       if (c1 && gm1->reuse() )               // must have type test first
  290.       { REPORT AddDS(gm1,gm2); gm2->tDelete(); return gm1; }
  291.       else if (c2 && gm2->reuse() )
  292.       { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
  293.       else
  294.       {
  295.          REPORT
  296.      GeneralMatrix* gmx = mtx.New(nr,nc,am); AddDS(gmx,gm1,gm2);
  297.      if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  298.          gmx->ReleaseAndDelete(); return gmx;
  299.       }
  300.    }
  301. }
  302.  
  303.  
  304. static GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
  305.    SubtractedMatrix* sm, MatrixType mtx)
  306. {
  307.    Tracer tr("GeneralSub");
  308.    if (!mtx) mtx = gm1->Type() + gm2->Type();
  309.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  310.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  311.       Throw(IncompatibleDimensionsException());
  312. #ifdef GXX
  313.    int c1,c2; sm->SelectVersion(mtx,c1,c2);
  314. #else
  315.    Boolean c1,c2; sm->SelectVersion(mtx,c1,c2); // causes problems for g++
  316. #endif
  317.    if (c1 && c2)
  318.    {
  319.       if (gm1->reuse())
  320.       { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
  321.       else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
  322.       else
  323.       {
  324.          REPORT
  325.      GeneralMatrix* gmx = gm1->Type().New(nr,nc,sm);
  326.          gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
  327.       }
  328.    }
  329.    else
  330.    {
  331.       if ( c1 && gm1->reuse() )
  332.       { REPORT  SubtractDS(gm1,gm2); gm2->tDelete(); return gm1; }
  333.       else if ( c2 && gm2->reuse() )
  334.       {
  335.          REPORT
  336.          ReverseSubtractDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
  337.       }
  338.       else
  339.       {
  340.          REPORT
  341.      GeneralMatrix* gmx = mtx.New(nr,nc,sm); SubtractDS(gmx,gm1,gm2);
  342.      if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  343.      gmx->ReleaseAndDelete(); return gmx;
  344.       }
  345.    }
  346. }
  347.  
  348. static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
  349.    MultipliedMatrix* mm, MatrixType mtx)
  350. {
  351.    REPORT
  352.    Tracer tr("GeneralMult1");
  353.    if (!mtx) mtx = gm1->Type() * gm2->Type();
  354.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  355.    if (gm1->Ncols() !=gm2->Nrows())
  356.       Throw(IncompatibleDimensionsException());
  357.    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  358.  
  359.    MatrixCol mcx(gmx, StoreOnExit+DirectPart);
  360.    MatrixCol mc2(gm2, LoadOnEntry);
  361.    while (nc--)
  362.    {
  363.       MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
  364.       Real* el = mcx();                              // pointer to an element
  365.       int n = mcx.Storage();
  366.       while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
  367.       mc2.Next(); mcx.Next();
  368.    }
  369.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  370. }
  371.  
  372. static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
  373.    MultipliedMatrix* mm, MatrixType mtx)
  374. {
  375.    REPORT
  376.    Tracer tr("GeneralMult2");
  377.    if (!mtx) mtx = gm1->Type() * gm2->Type();
  378.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  379.    if (gm1->Ncols() !=gm2->Nrows())
  380.       Throw(IncompatibleDimensionsException());
  381.    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  382.  
  383.    Real* el = gmx->Store(); int n = gmx->Storage();
  384.    while (n--) *el++ = 0.0;
  385.    MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
  386.    MatrixRow mr1(gm1, LoadOnEntry);
  387.    while (nr--)
  388.    {
  389.       MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
  390.       el = mr1();                              // pointer to an element
  391.       n = mr1.Storage();
  392.       while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
  393.       mr1.Next(); mrx.Next();
  394.    }
  395.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  396. }
  397.  
  398. static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
  399.    MultipliedMatrix* mm, MatrixType mtx)
  400. {
  401.    return GeneralMult1(gm1, gm2, mm, mtx);
  402. //   return GeneralMult2(gm1, gm2, mm, mtx);
  403. }
  404.  
  405. static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
  406.    BaseMatrix* sm, MatrixType mtx)
  407. {
  408.    REPORT
  409.    Tracer tr("GeneralSolv");
  410.    if (!mtx) mtx = gm1->Type().i()*gm2->Type();
  411.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  412.    if (gm1->Ncols() !=gm2->Nrows())
  413.       Throw(IncompatibleDimensionsException());
  414.    GeneralMatrix* gmx = mtx.New(nr,nc,sm);
  415.  
  416.    Real* r = new Real [nr]; MatrixErrorNoSpace(r);
  417. #ifndef ATandT
  418.    MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
  419.                            // deleted for ATandT, to balance deletion below
  420. #endif
  421.    MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
  422.    MatrixCol mc2(gm2, r, LoadOnEntry);
  423.    GeneralMatrix* gms = gm1->MakeSolver();
  424.    Try
  425.    {
  426.       while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
  427.    }
  428.    CatchAll
  429.    {
  430.       gms->tDelete(); delete gmx; gm2->tDelete();
  431. #ifndef ATandT
  432.       MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  433.                           // ATandT version 2.1 gives an internal error
  434. #endif
  435. #ifdef Version21
  436.       delete [] r;
  437. #else
  438.       delete [nr] r;
  439. #endif
  440.       Throw();
  441.    }
  442.    gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
  443. #ifndef ATandT
  444.    MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  445.                           // ATandT version 2.1 gives an internal error
  446. #endif
  447. #ifdef Version21
  448.    delete [] r;
  449. #else
  450.    delete [nr] r;
  451. #endif
  452.    return gmx;
  453. }
  454.  
  455. GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
  456. {
  457.    // Matrix Inversion - use solve routines
  458.    REPORT
  459.    gm=((BaseMatrix*&)bm)->Evaluate();
  460.    int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
  461. #ifdef TEMPS_DESTROYED_QUICKLY
  462.    GeneralMatrix* gmx;
  463.    Try { gmx = GeneralSolv(gm,&I,this,mtx); }
  464.    CatchAll { delete this; Throw(); }
  465.    delete this; return gmx;
  466. #else
  467.    return GeneralSolv(gm,&I,this,mtx);
  468. #endif   
  469. }
  470.  
  471.  
  472. /*************************** norm functions ****************************/
  473.  
  474. Real BaseMatrix::Norm1() const
  475. {
  476.    // maximum of sum of absolute values of a column
  477.    REPORT
  478.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  479.    int nc = gm->Ncols(); Real value = 0.0;
  480.    MatrixCol mc(gm, LoadOnEntry);
  481.    while (nc--)
  482.       { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
  483.    gm->tDelete(); return value;
  484. }
  485.  
  486. Real BaseMatrix::NormInfinity() const
  487. {
  488.    // maximum of sum of absolute values of a row
  489.    REPORT
  490.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  491.    int nr = gm->Nrows(); Real value = 0.0;
  492.    MatrixRow mr(gm, LoadOnEntry);
  493.    while (nr--)
  494.       { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
  495.    gm->tDelete(); return value;
  496. }
  497.