home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / NEWMAT8.CXX < prev    next >
C/C++ Source or Header  |  1992-11-29  |  9KB  |  342 lines
  1. //$$ newmat8.cxx         Advanced LU transform, scalar functions
  2.  
  3. // Copyright (C) 1991,2: R B Davies
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmatap.h"
  10.  
  11. //#define REPORT { static ExeCounter ExeCount(__LINE__,8); ExeCount++; }
  12.  
  13. #define REPORT {}
  14.  
  15.  
  16. /************************** LU transformation ****************************/
  17.  
  18. void CroutMatrix::ludcmp()
  19. // LU decomposition - from numerical recipes in C
  20. {
  21.    REPORT
  22.    Tracer trace("Crout(ludcmp)");
  23.    int i,j;
  24.  
  25.    Real* vv=new Real [nrows]; MatrixErrorNoSpace(vv);
  26.    MONITOR_REAL_NEW("Make  (CroutMat)",nrows,vv)
  27.    Real* a;
  28.  
  29.    a=store;
  30.    for (i=0;i<nrows;i++)
  31.    {
  32.       Real big=0.0;
  33.       j=nrows; while (j--) { Real temp=fabs(*a++); if (temp > big) big=temp; }
  34.       if (big == 0.0)
  35.       {
  36.          MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  37. #ifdef Version21
  38.          sing = TRUE; delete [] vv; return;
  39. #else
  40.          sing = TRUE; delete [nrows] vv; return;
  41. #endif
  42.       }
  43.       vv[i]=1.0/big;
  44.    }
  45.  
  46.    Real* aj=store;
  47.    for (j=0;j<nrows;j++)
  48.    {
  49.       Real* ai=store;
  50.       for (i=0;i<j;i++)
  51.       {
  52.          Real sum=*(ai+j); Real* aix=ai; Real* ajx=aj;
  53.          int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  54.          *ajx = sum; ai += nrows;
  55.       }
  56.  
  57.       Real big=0.0; int imax;
  58.       for (i=j;i<nrows;i++)
  59.       {
  60.          Real sum=*(ai+j); Real* aix=ai; Real* ajx=aj;
  61.          int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  62.          *aix = sum; ai += nrows;
  63.          Real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
  64.       }
  65.  
  66.       if (j != imax)
  67.       {
  68.          Real* amax=store+imax*nrows; Real* ajx=store+j*nrows;
  69.          int k=nrows; while(k--) { Real dum=*amax; *amax++=*ajx; *ajx++=dum; }
  70.          d=!d; vv[imax]=vv[j];
  71.       }
  72.  
  73.       indx[j]=imax; ai=aj+j*nrows;
  74.       if (*ai == 0.0)
  75.       {
  76.          MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  77. #ifdef Version21
  78.          sing = TRUE; delete [] vv; return;
  79. #else
  80.          sing = TRUE; delete [nrows] vv; return;
  81. #endif
  82.       }
  83.       Real dum=1.0/(*ai);
  84.       i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
  85.  
  86.       aj++;
  87.    }
  88.    MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  89. #ifdef Version21
  90.    delete [] vv;
  91. #else
  92.    delete [nrows] vv;
  93. #endif
  94. }
  95.  
  96. void CroutMatrix::lubksb(Real* B, int mini)
  97. {
  98.    REPORT
  99.    Tracer trace("Crout(lubksb)");
  100.    if (sing) Throw(SingularException(*this));   
  101.    int i,j; int ii=-1; Real* ai=store;
  102.  
  103.    for (i=0;i<nrows;i++)
  104.    {
  105.       int ip=indx[i]; Real sum=B[ip]; B[ip]=B[i];
  106.       if (ii>=0)
  107.       {
  108.          Real* aij=ai+ii; Real* bj=B+ii; j=i-ii;
  109.          while (j--) { sum -= *aij++ * *bj; bj++; }
  110.       }
  111.       else if (sum) ii=i;
  112.       B[i]=sum; ai += nrows;
  113.    }
  114.  
  115.    for (i=nrows-1;i>=mini;i--)
  116.    {
  117.       Real* bj=B+i; ai -= nrows; Real* ajx=ai+i; Real sum=*bj; bj++;
  118.       j=nrows-i; while(--j) { sum -= *(++ajx) * *bj; bj++; }
  119.       B[i]=sum/(*(ai+i));
  120.    }
  121. }
  122.  
  123.  
  124. /****************************** scalar functions ****************************/
  125.  
  126. inline Real square(Real x) { return x*x; }
  127.  
  128. Real GeneralMatrix::SumSquare() const
  129. {
  130.    REPORT
  131.    Real sum = 0.0; int i = storage; Real* s = store;
  132.    while (i--) sum += square(*s++);
  133.    ((GeneralMatrix&)*this).tDelete(); return sum;
  134. }
  135.  
  136. Real GeneralMatrix::SumAbsoluteValue() const
  137. {
  138.    REPORT
  139.    Real sum = 0.0; int i = storage; Real* s = store;
  140.    while (i--) sum += fabs(*s++);
  141.    ((GeneralMatrix&)*this).tDelete(); return sum;
  142. }
  143.  
  144. Real GeneralMatrix::MaximumAbsoluteValue() const
  145. {
  146.    REPORT
  147.    Real maxval = 0.0; int i = storage; Real* s = store;
  148.    while (i--) { Real a = fabs(*s++); if (maxval < a) maxval = a; }
  149.    ((GeneralMatrix&)*this).tDelete(); return maxval;
  150. }
  151.  
  152. Real SymmetricMatrix::SumSquare() const
  153. {
  154.    REPORT
  155.    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
  156.    for (int i = 0; i<nr; i++)
  157.    {
  158.       int j = i;
  159.       while (j--) sum2 += square(*s++);
  160.       sum1 += square(*s++);
  161.    }
  162.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  163. }
  164.  
  165. Real SymmetricMatrix::SumAbsoluteValue() const
  166. {
  167.    REPORT
  168.    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
  169.    for (int i = 0; i<nr; i++)
  170.    {
  171.       int j = i;
  172.       while (j--) sum2 += fabs(*s++);
  173.       sum1 += fabs(*s++);
  174.    }
  175.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  176. }
  177.  
  178. Real BaseMatrix::SumSquare() const
  179. {
  180.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  181.    Real s = gm->SumSquare(); return s;
  182. }
  183.  
  184. Real BaseMatrix::SumAbsoluteValue() const
  185. {
  186.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  187.    Real s = gm->SumAbsoluteValue(); return s;
  188. }
  189.  
  190. Real BaseMatrix::MaximumAbsoluteValue() const
  191. {
  192.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  193.    Real s = gm->MaximumAbsoluteValue(); return s;
  194. }
  195.  
  196. Real Matrix::Trace() const
  197. {
  198.    REPORT
  199.    Tracer trace("Trace");
  200.    int i = nrows; int d = i+1;
  201.    if (i != ncols) Throw(NotSquareException(*this));
  202.    Real sum = 0.0; Real* s = store;
  203.    while (i--) { sum += *s; s += d; }
  204.    ((GeneralMatrix&)*this).tDelete(); return sum;
  205. }
  206.  
  207. Real DiagonalMatrix::Trace() const
  208. {
  209.    REPORT
  210.    int i = nrows; Real sum = 0.0; Real* s = store;
  211.    while (i--) sum += *s++;
  212.    ((GeneralMatrix&)*this).tDelete(); return sum;
  213. }
  214.  
  215. Real SymmetricMatrix::Trace() const
  216. {
  217.    REPORT
  218.    int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
  219.    while (i--) { sum += *s; s += j++; }
  220.    ((GeneralMatrix&)*this).tDelete(); return sum;
  221. }
  222.  
  223. Real LowerTriangularMatrix::Trace() const
  224. {
  225.    REPORT
  226.    int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
  227.    while (i--) { sum += *s; s += j++; }
  228.    ((GeneralMatrix&)*this).tDelete(); return sum;
  229. }
  230.  
  231. Real UpperTriangularMatrix::Trace() const
  232. {
  233.    REPORT
  234.    int i = nrows; Real sum = 0.0; Real* s = store;
  235.    while (i) { sum += *s; s += i--; }
  236.    ((GeneralMatrix&)*this).tDelete(); return sum;
  237. }
  238.  
  239. Real BandMatrix::Trace() const
  240. {
  241.    REPORT
  242.    int i = nrows; int w = lower+upper+1;
  243.    Real sum = 0.0; Real* s = store+lower;
  244.    while (i--) { sum += *s; s += w; }
  245.    ((GeneralMatrix&)*this).tDelete(); return sum;
  246. }
  247.  
  248. Real SymmetricBandMatrix::Trace() const
  249. {
  250.    REPORT
  251.    int i = nrows; int w = lower+1;
  252.    Real sum = 0.0; Real* s = store+lower;
  253.    while (i--) { sum += *s; s += w; }
  254.    ((GeneralMatrix&)*this).tDelete(); return sum;
  255. }
  256.  
  257. Real BaseMatrix::Trace() const
  258. {
  259.    REPORT
  260.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(MatrixType::Dg);
  261.    Real sum = gm->Trace(); return sum;
  262. }
  263.  
  264. void LogAndSign::operator*=(Real x)
  265. {
  266.    if (x > 0.0) { log_value += log(x); }
  267.    else if (x < 0.0) { log_value += log(-x); sign = -sign; }
  268.    else sign = 0;
  269. }
  270.  
  271. Real LogAndSign::Value() const { return sign * exp(log_value); }
  272.  
  273. LogAndSign::LogAndSign(Real f)
  274. {
  275.    if (f == 0.0) { log_value = 0.0; sign = 0; return; }
  276.    else if (f < 0.0) { sign = -1; f = -f; }
  277.    else sign = 1;
  278.    log_value = log(f);
  279. }
  280.  
  281. LogAndSign DiagonalMatrix::LogDeterminant() const
  282. {
  283.    REPORT
  284.    int i = nrows; LogAndSign sum; Real* s = store;
  285.    while (i--) sum *= *s++;
  286.    ((GeneralMatrix&)*this).tDelete(); return sum;
  287. }
  288.  
  289. LogAndSign LowerTriangularMatrix::LogDeterminant() const
  290. {
  291.    REPORT
  292.    int i = nrows; LogAndSign sum; Real* s = store; int j = 2;
  293.    while (i--) { sum *= *s; s += j++; }
  294.    ((GeneralMatrix&)*this).tDelete(); return sum;
  295. }
  296.  
  297. LogAndSign UpperTriangularMatrix::LogDeterminant() const
  298. {
  299.    REPORT
  300.    int i = nrows; LogAndSign sum; Real* s = store;
  301.    while (i) { sum *= *s; s += i--; }
  302.    ((GeneralMatrix&)*this).tDelete(); return sum;
  303. }
  304.  
  305. LogAndSign BaseMatrix::LogDeterminant() const
  306. {
  307.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  308.    LogAndSign sum = gm->LogDeterminant(); return sum;
  309. }
  310.  
  311. LogAndSign GeneralMatrix::LogDeterminant() const
  312. {
  313.    REPORT
  314.    Tracer tr("Determinant");
  315.    if (nrows != ncols) Throw(NotSquareException(*this));
  316.    CroutMatrix C(*this); return C.LogDeterminant();
  317. }
  318.  
  319. LogAndSign CroutMatrix::LogDeterminant() const
  320. {
  321.    REPORT
  322.    if (sing) return 0.0;
  323.    int i = nrows; int dd = i+1; LogAndSign sum; Real* s = store;
  324.    while (i--) { sum *= *s; s += dd; }
  325.    if (!d) sum.ChangeSign(); return sum;
  326.  
  327. }
  328.  
  329. LinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm)
  330. : gm( ( ((BaseMatrix&)bm).Evaluate() )->MakeSolver() )
  331. {
  332.    if (gm==&bm)
  333.    // want a copy if  *gm is actually bm
  334.    {
  335.       REPORT
  336.       GeneralMatrix* gmx = gm->Type().New();
  337.       gmx->GetMatrix(gm); gm = gmx;
  338.    } 
  339.    else { REPORT  gm->Protect(); }
  340. }
  341.  
  342.