home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / newmat03 / newmat7.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-12  |  11.6 KB  |  366 lines
  1. //$$ newmat7.cxx     Invert, solve, binary operations
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "newmat.hxx"
  8. #include "newmatrc.hxx"
  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 void Add(GeneralMatrix*,GeneralMatrix*, GeneralMatrix*);
  83. static void Add(GeneralMatrix*,GeneralMatrix*);
  84.  
  85. static void Subtract(GeneralMatrix*,GeneralMatrix*, GeneralMatrix*);
  86. static void Subtract(GeneralMatrix*,GeneralMatrix*);
  87. static void ReverseSubtract(GeneralMatrix*,GeneralMatrix*);
  88.  
  89. static GeneralMatrix* GeneralAdd(GeneralMatrix*,GeneralMatrix*,MatrixType);
  90. static GeneralMatrix* GeneralSub(GeneralMatrix*,GeneralMatrix*,MatrixType);
  91. static GeneralMatrix* GeneralMult(GeneralMatrix*,GeneralMatrix*,MatrixType);
  92. static GeneralMatrix* GeneralSolv(GeneralMatrix*,GeneralMatrix*,MatrixType);
  93.  
  94. GeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
  95. { REPORT return GeneralAdd(bm1->Evaluate(), bm2->Evaluate(), mt); }
  96.  
  97. GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
  98. { REPORT return GeneralSub(bm1->Evaluate(), bm2->Evaluate(), mt); }
  99.  
  100. GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
  101. {
  102.    REPORT
  103.    GeneralMatrix* gm2 = bm2->Evaluate();
  104.    MatrixType mtsym(MatrixType::Sym);
  105.    if (gm2->Type() == mtsym) gm2 = gm2->Evaluate(MatrixType::Rect);
  106.    return GeneralMult(bm1->Evaluate(), gm2, mt);
  107. }
  108.  
  109. GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
  110. { REPORT return GeneralSolv(bm1->Evaluate(), bm2->Evaluate(), mt); }
  111.  
  112. // routines for adding or subtracting matrices of identical storage structure
  113.  
  114. static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  115. {
  116.    REPORT
  117.    register real* s1=gm1->Store(); register real* s2=gm2->Store();
  118.    register real* s=gm->Store(); register int i=gm->Storage();
  119.    while (i--) *s++ = *s1++ + *s2++;
  120. }
  121.    
  122. static void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
  123. {
  124.    REPORT
  125.    register real* s2=gm2->Store();
  126.    register real* s=gm->Store(); register int i=gm->Storage();
  127.    while (i--) *s++ += *s2++;
  128. }
  129.  
  130. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  131. {
  132.    REPORT
  133.    register real* s1=gm1->Store(); register real* s2=gm2->Store();
  134.    register real* s=gm->Store(); register int i=gm->Storage();
  135.    while (i--) *s++ = *s1++ - *s2++;
  136. }
  137.  
  138. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  139. {
  140.    REPORT
  141.    register real* s2=gm2->Store();
  142.    register real* s=gm->Store(); register int i=gm->Storage();
  143.    while (i--) *s++ -= *s2++;
  144. }
  145.  
  146. static void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  147. {
  148.    REPORT
  149.    register real* s2=gm2->Store();
  150.    register real* s=gm->Store(); register int i=gm->Storage();
  151.    while (i--) { *s = *s2++ - *s; s++; }
  152. }
  153.  
  154. static GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
  155.    MatrixType mtx)
  156. {
  157.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  158.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  159.       MatrixError("Incompatible dimensions");
  160.    if (!mtx) mtx = gm1->Type() + gm2->Type();
  161.    if (mtx == gm1->Type() && mtx == gm2->Type())
  162.                          // modify when band or sparse
  163.                          // matrices are included.
  164.    {
  165.       if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
  166.       else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
  167.       else
  168.       {
  169.          REPORT
  170.          GeneralMatrix* gmx = gm1->Type().New(nr,nc); gmx->ReleaseAndDelete();
  171.          gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
  172.       }
  173.    }
  174.    else
  175.    {
  176.       if (gm1->Type()==mtx && gm1->reuse() )    // must have type test first
  177.       {
  178.      REPORT
  179.      MatrixRow mr1(gm1, StoreOnExit+LoadOnEntry+DirectPart);
  180.      MatrixRow mr2(gm2, LoadOnEntry);
  181.      while (nr--) { mr1.Add(mr2); mr1.Next(); mr2.Next(); }
  182.          gm2->tDelete(); return gm1;
  183.       }
  184.       else if (gm2->Type()==mtx && gm2->reuse() )
  185.       {
  186.          REPORT
  187.      MatrixRow mr1(gm1, LoadOnEntry);
  188.      MatrixRow mr2(gm2, StoreOnExit+LoadOnEntry+DirectPart);
  189.      while (nr--) { mr2.Add(mr1); mr1.Next(); mr2.Next(); }
  190.      if (gm1->Type()!=mtx) gm1->tDelete();
  191.          return gm2;
  192.       }
  193.       else
  194.       {
  195.          REPORT
  196.      GeneralMatrix* gmx = mtx.New(nr,nc);
  197.          MatrixRow mr1(gm1, LoadOnEntry);
  198.      MatrixRow mr2(gm2, LoadOnEntry);
  199.      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  200.      while (nr--)
  201.      { mrx.Add(mr1,mr2); mr1.Next(); mr2.Next(); mrx.Next(); }
  202.      if (gm1->Type()!=mtx) gm1->tDelete();
  203.      if (gm2->Type()!=mtx) gm2->tDelete();
  204.          gmx->ReleaseAndDelete(); return gmx;
  205.       }
  206.    }
  207. }
  208.  
  209.  
  210. static GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
  211.    MatrixType mtx)
  212. {
  213.    if (!mtx) mtx = gm1->Type() - gm2->Type();
  214.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  215.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  216.       MatrixError("Incompatible dimensions");
  217.    if (mtx == gm1->Type() && mtx == gm2->Type())
  218.                          // modify when band or sparse
  219.                                              // matrices are included.
  220.    {
  221.       if (gm1->reuse())
  222.       { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
  223.       else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
  224.       else
  225.       {
  226.          REPORT
  227.      GeneralMatrix* gmx = gm1->Type().New(nr,nc); gmx->ReleaseAndDelete();
  228.          gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
  229.       }
  230.    }
  231.    else
  232.    {
  233.       if (gm1->Type() == mtx && gm1->reuse() )
  234.       {
  235.          REPORT
  236.      MatrixRow mr1(gm1, LoadOnEntry+StoreOnExit+DirectPart);
  237.      MatrixRow mr2(gm2, LoadOnEntry);
  238.      while (nr--) { mr1.Sub(mr2); mr1.Next(); mr2.Next(); }
  239.          gm2->tDelete(); return gm1;
  240.       }
  241.       else if (gm2->Type() == mtx && gm2->reuse() )
  242.       {
  243.          REPORT
  244.          MatrixRow mr1(gm1, LoadOnEntry);
  245.      MatrixRow mr2(gm2, LoadOnEntry+StoreOnExit+DirectPart);
  246.      while (nr--) { mr2.RevSub(mr1); mr1.Next(); mr2.Next(); }
  247.      if (gm1->Type()!=mtx) gm1->tDelete();
  248.      return gm2;
  249.       }
  250.       else
  251.       {
  252.          REPORT
  253.      GeneralMatrix* gmx = mtx.New(nr,nc);
  254.          MatrixRow mr1(gm1, LoadOnEntry);
  255.          MatrixRow mr2(gm2, LoadOnEntry);
  256.      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  257.      while (nr--)
  258.         { mrx.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mrx.Next(); }
  259.      if (gm1->Type()!=mtx) gm1->tDelete();
  260.      if (gm2->Type()!=mtx) gm2->tDelete();
  261.      gmx->ReleaseAndDelete(); return gmx;
  262.       }
  263.    }
  264. }
  265.  
  266. /*
  267. static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
  268.    MatrixType mtx)
  269. {
  270.    REPORT
  271.    if (!mtx) mtx = gm1->Type() * gm2->Type();
  272.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  273.    if (gm1->Ncols() !=gm2->Nrows()) MatrixError("Incompatible dimensions");
  274.    GeneralMatrix* gmx = mtx.New(nr,nc);
  275.  
  276.    MatrixCol mcx(gmx, StoreOnExit+DirectPart);
  277.    MatrixCol mc2(gm2, LoadOnEntry);
  278.    while (nc--)
  279.    {
  280.       MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
  281.       real* el = mcx();                              // pointer to an element
  282.       int n = mcx.Storage();
  283.       while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
  284.       mc2.Next(); mcx.Next();
  285.    }
  286.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  287. }
  288. */
  289.  
  290. static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
  291.    MatrixType mtx)
  292. {
  293.    REPORT
  294.    if (!mtx) mtx = gm1->Type() * gm2->Type();
  295.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  296.    if (gm1->Ncols() !=gm2->Nrows()) MatrixError("Incompatible dimensions");
  297.    GeneralMatrix* gmx = mtx.New(nr,nc);
  298.  
  299.    real* el = gmx->Store(); int n = gmx->Storage();
  300.    while (n--) *el++ = 0.0;
  301.    MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
  302.    MatrixRow mr1(gm1, LoadOnEntry);
  303.    while (nr--)
  304.    {
  305.       MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
  306.       el = mr1();                              // pointer to an element
  307.       n = mr1.Storage();
  308.       while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
  309.       mr1.Next(); mrx.Next();
  310.    }
  311.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  312. }
  313.  
  314. static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
  315.    MatrixType mtx)
  316. {
  317.    REPORT
  318.    if (!mtx) mtx = gm1->Type().i()*gm2->Type();
  319.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  320.    if (gm1->Ncols() !=gm2->Nrows()) MatrixError("Incompatible dimensions");
  321.    GeneralMatrix* gmx = mtx.New(nr,nc);
  322.  
  323.    real* r = new real [nr]; MatrixErrorNoSpace(r);
  324.    MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
  325.    MatrixCol mc2(gm2, r, LoadOnEntry);
  326.    GeneralMatrix* gms = gm1->MakeSolver();
  327.    while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
  328.    gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
  329.    delete [nr] r;
  330.    return gmx;
  331. }
  332.  
  333. GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
  334. {
  335.    // Matrix Inversion - use solve routines
  336.    REPORT
  337.    GeneralMatrix* gm=bm->Evaluate();
  338.    int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
  339.    return GeneralSolv(gm,&I,mtx);
  340. }
  341.  
  342.  
  343. /*************************** norm functions ****************************/
  344.  
  345. real BaseMatrix::Norm1()
  346. {
  347.    // maximum of sum of absolute values of a column
  348.    REPORT
  349.    GeneralMatrix* gm = Evaluate(); int nc = gm->Ncols(); real value = 0.0;
  350.    MatrixCol mc(gm, LoadOnEntry);
  351.    while (nc--)
  352.       { real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
  353.    gm->tDelete(); return value;
  354. }
  355.  
  356. real BaseMatrix::NormInfinity()
  357. {
  358.    // maximum of sum of absolute values of a row
  359.    REPORT
  360.    GeneralMatrix* gm = Evaluate(); int nr = gm->Nrows(); real value = 0.0;
  361.    MatrixRow mr(gm, LoadOnEntry);
  362.    while (nr--)
  363.       { real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
  364.    gm->tDelete(); return value;
  365. }
  366.