home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / sources / misc / 4253 / newmat7.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  15.7 KB  |  540 lines
  1. //$$ newmat7.cxx     Invert, solve, binary operations
  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. //#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;
  98.    Try { gmx = GeneralAdd(gm1,gm2,this,mt); }
  99.    CatchAll { delete this; Throw(); }
  100.    delete this; return gmx;
  101. #else
  102.    return GeneralAdd(gm1,gm2,this,mt);
  103. #endif   
  104. }
  105.  
  106. GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
  107. {
  108.    REPORT
  109.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  110.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  111. #ifdef TEMPS_DESTROYED_QUICKLY
  112.    GeneralMatrix* gmx;
  113.    Try { gmx = GeneralSub(gm1,gm2,this,mt); }
  114.    CatchAll { delete this; Throw(); }
  115.    delete this; return gmx;
  116. #else
  117.    return GeneralSub(gm1,gm2,this,mt);
  118. #endif   
  119. }
  120.  
  121. GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
  122. {
  123.    REPORT
  124.    gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  125.    gm2 = gm2->Evaluate(gm2->Type().MultRHS());     // no symmetric on RHS
  126.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  127. #ifdef TEMPS_DESTROYED_QUICKLY
  128.    GeneralMatrix* gmx;
  129.    Try { gmx = GeneralMult(gm1, gm2, this, mt); }
  130.    CatchAll { delete this; Throw(); }
  131.    delete this; return gmx;
  132. #else
  133.    return GeneralMult(gm1, gm2, this, mt);
  134. #endif   
  135. }
  136.  
  137. GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
  138. {
  139.    REPORT
  140.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  141.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  142. #ifdef TEMPS_DESTROYED_QUICKLY
  143.    GeneralMatrix* gmx;
  144.    Try { gmx = GeneralSolv(gm1,gm2,this,mt); }
  145.    CatchAll { delete this; Throw(); }
  146.    delete this; return gmx;
  147. #else
  148.    return GeneralSolv(gm1,gm2,this,mt);
  149. #endif   
  150. }
  151.  
  152. // routines for adding or subtracting matrices of identical storage structure
  153.  
  154. static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  155. {
  156.    REPORT
  157.    Real* s1=gm1->Store(); Real* s2=gm2->Store();
  158.    Real* s=gm->Store(); int i=gm->Storage() >> 2;
  159.    while (i--)
  160.    {
  161.        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  162.        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  163.    }
  164.    i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
  165. }
  166.    
  167. static void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
  168. {
  169.    REPORT
  170.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  171.    while (i--)
  172.    { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
  173.    i=gm->Storage() & 3; while (i--) *s++ += *s2++;
  174. }
  175.  
  176. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  177. {
  178.    REPORT
  179.    Real* s1=gm1->Store(); Real* s2=gm2->Store();
  180.    Real* s=gm->Store(); int i=gm->Storage() >> 2;
  181.    while (i--)
  182.    {
  183.        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  184.        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  185.    }
  186.    i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
  187. }
  188.  
  189. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  190. {
  191.    REPORT
  192.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  193.    while (i--)
  194.    { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
  195.    i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
  196. }
  197.  
  198. static void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  199. {
  200.    REPORT
  201.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  202.    while (i--)
  203.    {
  204.       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  205.       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  206.    }
  207.    i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
  208. }
  209.  
  210. // routines for adding or subtracting matrices of different storage structure
  211.  
  212. static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  213. {
  214.    int nr = gm->Nrows();
  215.    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  216.    MatrixRow mr(gm, StoreOnExit+DirectPart);
  217.    while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  218. }
  219.    
  220. static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  221. // Add into first argument
  222. {
  223.    int nr = gm->Nrows();
  224.    MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
  225.    MatrixRow mr2(gm2, LoadOnEntry);
  226.    while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
  227. }
  228.  
  229. static void SubtractDS
  230.    (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  231. {
  232.    int nr = gm->Nrows();
  233.    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  234.    MatrixRow mr(gm, StoreOnExit+DirectPart);
  235.    while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  236. }
  237.  
  238. static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  239. {
  240.    int nr = gm->Nrows();
  241.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  242.    MatrixRow mr2(gm2, LoadOnEntry);
  243.    while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
  244. }
  245.  
  246. static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  247. {
  248.    int nr = gm->Nrows();
  249.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  250.    MatrixRow mr2(gm2, LoadOnEntry);
  251.    while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
  252. }
  253.  
  254. #ifdef __GNUG__
  255. void AddedMatrix::SelectVersion
  256.    (MatrixType mtx, int& c1, int& c2) const
  257. #else
  258. void AddedMatrix::SelectVersion
  259.    (MatrixType mtx, Boolean& c1, Boolean& c2) const
  260. #endif
  261. // for determining version of add and subtract routines
  262. // will need to modify if further matrix structures are introduced
  263. {
  264.    MatrixBandWidth bm1 = gm1->BandWidth();
  265.    MatrixBandWidth bm2 = gm2->BandWidth();
  266.    MatrixBandWidth bmx = bm1 + bm2;
  267.    c1 = (mtx == gm1->Type()) && (bmx == bm1);
  268.    c2 = (mtx == gm2->Type()) && (bmx == bm2);
  269. }
  270.  
  271. static GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
  272.    AddedMatrix* am, MatrixType mtx)
  273. {
  274.    Tracer tr("GeneralAdd");
  275.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  276.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols()) 
  277.       Throw(IncompatibleDimensionsException());
  278.    Compare(gm1->Type() + gm2->Type(),mtx);
  279. #ifdef __GNUG__
  280.    int c1,c2; am->SelectVersion(mtx,c1,c2);
  281. #else
  282.    Boolean c1,c2; am->SelectVersion(mtx,c1,c2); // causes problems for g++
  283. #endif
  284.    if (c1 && c2)
  285.    {
  286.       if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
  287.       else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
  288.       else
  289.       {
  290.          REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
  291.          gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
  292.       }
  293.    }
  294.    else
  295.    {
  296.       if (c1 && gm1->reuse() )               // must have type test first
  297.       { REPORT AddDS(gm1,gm2); gm2->tDelete(); return gm1; }
  298.       else if (c2 && gm2->reuse() )
  299.       { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
  300.       else
  301.       {
  302.          REPORT
  303.      GeneralMatrix* gmx = mtx.New(nr,nc,am); AddDS(gmx,gm1,gm2);
  304.      if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  305.          gmx->ReleaseAndDelete(); return gmx;
  306.       }
  307.    }
  308. }
  309.  
  310.  
  311. static GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
  312.    SubtractedMatrix* sm, MatrixType mtx)
  313. {
  314.    Tracer tr("GeneralSub");
  315.    Compare(gm1->Type() + gm2->Type(),mtx);
  316.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  317.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  318.       Throw(IncompatibleDimensionsException());
  319. #ifdef __GNUG__
  320.    int c1,c2; sm->SelectVersion(mtx,c1,c2);
  321. #else
  322.    Boolean c1,c2; sm->SelectVersion(mtx,c1,c2); // causes problems for g++
  323. #endif
  324.    if (c1 && c2)
  325.    {
  326.       if (gm1->reuse())
  327.       { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
  328.       else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
  329.       else
  330.       {
  331.          REPORT
  332.      GeneralMatrix* gmx = gm1->Type().New(nr,nc,sm);
  333.          gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
  334.       }
  335.    }
  336.    else
  337.    {
  338.       if ( c1 && gm1->reuse() )
  339.       { REPORT  SubtractDS(gm1,gm2); gm2->tDelete(); return gm1; }
  340.       else if ( c2 && gm2->reuse() )
  341.       {
  342.          REPORT
  343.          ReverseSubtractDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
  344.       }
  345.       else
  346.       {
  347.          REPORT
  348.      GeneralMatrix* gmx = mtx.New(nr,nc,sm); SubtractDS(gmx,gm1,gm2);
  349.      if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  350.      gmx->ReleaseAndDelete(); return gmx;
  351.       }
  352.    }
  353. }
  354.  
  355. static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
  356.    MultipliedMatrix* mm, MatrixType mtx)
  357. {
  358.    REPORT
  359.    Tracer tr("GeneralMult1");
  360.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  361.    if (gm1->Ncols() !=gm2->Nrows())
  362.       Throw(IncompatibleDimensionsException());
  363.    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  364.  
  365.    MatrixCol mcx(gmx, StoreOnExit+DirectPart);
  366.    MatrixCol mc2(gm2, LoadOnEntry);
  367.    while (nc--)
  368.    {
  369.       MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
  370.       Real* el = mcx();                              // pointer to an element
  371.       int n = mcx.Storage();
  372.       while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
  373.       mc2.Next(); mcx.Next();
  374.    }
  375.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  376. }
  377.  
  378. static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
  379.    MultipliedMatrix* mm, MatrixType mtx)
  380. {
  381.    // version that accesses by row only - not good for thin matrices
  382.    // or column vectors in right hand term. Needs fixing
  383.    REPORT
  384.    Tracer tr("GeneralMult2");
  385.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  386.    if (gm1->Ncols() !=gm2->Nrows())
  387.       Throw(IncompatibleDimensionsException());
  388.    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  389.  
  390.    Real* el = gmx->Store(); int n = gmx->Storage();
  391.    while (n--) *el++ = 0.0;
  392.    MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
  393.    MatrixRow mr1(gm1, LoadOnEntry);
  394.    while (nr--)
  395.    {
  396.       MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
  397.       el = mr1();                              // pointer to an element
  398.       n = mr1.Storage();
  399.       while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
  400.       mr1.Next(); mrx.Next();
  401.    }
  402.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  403. }
  404.  
  405. static GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
  406. {
  407.    // matrix multiplication for type Matrix only
  408.    REPORT
  409.    Tracer tr("MatrixMult");
  410.  
  411.    int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
  412.    if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
  413.  
  414.    Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
  415.  
  416.    Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
  417.    
  418.    if (ncr)
  419.    {
  420.       while (nr--)
  421.       {
  422.          Real* s2x = s2; int j = ncr;
  423.          Real* sx = s; Real f = *s1++; int k = nc;
  424.          while (k--) *sx++ = f * *s2x++;
  425.          while (--j)
  426.             { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
  427.          s = sx;
  428.       }
  429.    }
  430.    else *gm = 0.0;
  431.  
  432.    gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
  433. }
  434.  
  435. static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
  436.    MultipliedMatrix* mm, MatrixType mtx)
  437. {
  438.    if ( Rectangular(gm1->Type(), gm2->Type(), mtx)) return mmMult(gm1, gm2);
  439.    else
  440.    {
  441.       Compare(gm1->Type() * gm2->Type(),mtx);
  442.       int nr = gm2->Nrows(); int nc = gm2->Ncols();
  443.       if (nc <= 5 && nr > nc) return GeneralMult1(gm1, gm2, mm, mtx);
  444.       else return GeneralMult2(gm1, gm2, mm, mtx);
  445.    }
  446. }
  447.  
  448. static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
  449.    BaseMatrix* sm, MatrixType mtx)
  450. {
  451.    REPORT
  452.    Tracer tr("GeneralSolv");
  453.    Compare(gm1->Type().i() * gm2->Type(),mtx);
  454.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  455.    if (gm1->Ncols() !=gm2->Nrows())
  456.       Throw(IncompatibleDimensionsException());
  457.    GeneralMatrix* gmx = mtx.New(nr,nc,sm);
  458.  
  459.    Real* r = new Real [nr]; MatrixErrorNoSpace(r);
  460. #ifndef ATandT
  461.    MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
  462.                            // deleted for ATandT, to balance deletion below
  463. #endif
  464.    MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
  465.    MatrixCol mc2(gm2, r, LoadOnEntry);
  466.    GeneralMatrix* gms = gm1->MakeSolver();
  467.    Try
  468.    {
  469.       while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
  470.    }
  471.    CatchAll
  472.    {
  473.       gms->tDelete(); delete gmx; gm2->tDelete();
  474. #ifndef ATandT
  475.       MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  476.                           // ATandT version 2.1 gives an internal error
  477. #endif
  478. #ifdef Version21
  479.       delete [] r;
  480. #else
  481.       delete [nr] r;
  482. #endif
  483.       Throw();
  484.    }
  485.    gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
  486. #ifndef ATandT
  487.    MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  488.                           // ATandT version 2.1 gives an internal error
  489. #endif
  490. #ifdef Version21
  491.    delete [] r;
  492. #else
  493.    delete [nr] r;
  494. #endif
  495.    return gmx;
  496. }
  497.  
  498. GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
  499. {
  500.    // Matrix Inversion - use solve routines
  501.    REPORT
  502.    gm=((BaseMatrix*&)bm)->Evaluate();
  503.    int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
  504. #ifdef TEMPS_DESTROYED_QUICKLY
  505.    GeneralMatrix* gmx;
  506.    Try { gmx = GeneralSolv(gm,&I,this,mtx); }
  507.    CatchAll { delete this; Throw(); }
  508.    delete this; return gmx;
  509. #else
  510.    return GeneralSolv(gm,&I,this,mtx);
  511. #endif   
  512. }
  513.  
  514.  
  515. /*************************** norm functions ****************************/
  516.  
  517. Real BaseMatrix::Norm1() const
  518. {
  519.    // maximum of sum of absolute values of a column
  520.    REPORT
  521.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  522.    int nc = gm->Ncols(); Real value = 0.0;
  523.    MatrixCol mc(gm, LoadOnEntry);
  524.    while (nc--)
  525.       { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
  526.    gm->tDelete(); return value;
  527. }
  528.  
  529. Real BaseMatrix::NormInfinity() const
  530. {
  531.    // maximum of sum of absolute values of a row
  532.    REPORT
  533.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  534.    int nr = gm->Nrows(); Real value = 0.0;
  535.    MatrixRow mr(gm, LoadOnEntry);
  536.    while (nr--)
  537.       { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
  538.    gm->tDelete(); return value;
  539. }
  540.