home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / newmat04 / part04 / newmat8.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-30  |  6.1 KB  |  274 lines
  1. //$$ newmat8.cxx         Advanced LU transform, scalar functions
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.hxx"
  8.  
  9. #include "newmatap.hxx"
  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.    int i,j;
  23.  
  24.    real* vv=new real [nrows]; MatrixErrorNoSpace(vv);
  25.    real* a;
  26.  
  27.    a=store;
  28.    for (i=0;i<nrows;i++)
  29.    {
  30.       real big=0.0;
  31.       j=nrows; while (j--) { real temp=fabs(*a++); if (temp > big) big=temp; }
  32.       if (big == 0.0) MatrixError("Singular matrix");
  33.       vv[i]=1.0/big;
  34.    }
  35.  
  36.    real* aj=store;
  37.    for (j=0;j<nrows;j++)
  38.    {
  39.       real* ai=store;
  40.       for (i=0;i<j;i++)
  41.       {
  42.          real sum=*(ai+j); real* aix=ai; real* ajx=aj;
  43.          int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  44.          *ajx = sum; ai += nrows;
  45.       }
  46.  
  47.       real big=0.0; int imax;
  48.       for (i=j;i<nrows;i++)
  49.       {
  50.          real sum=*(ai+j); real* aix=ai; real* ajx=aj;
  51.          int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  52.          *aix = sum; ai += nrows;
  53.          real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
  54.       }
  55.  
  56.       if (j != imax)
  57.       {
  58.          real* amax=store+imax*nrows; real* ajx=store+j*nrows;
  59.          int k=nrows; while(k--) { real dum=*amax; *amax++=*ajx; *ajx++=dum; }
  60.          d=!d; vv[imax]=vv[j];
  61.       }
  62.  
  63.       indx[j]=imax; ai=aj+j*nrows;
  64.       if (*ai == 0.0) MatrixError("Singular matrix");
  65.       real dum=1.0/(*ai);
  66.       i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
  67.  
  68.       aj++;
  69.    }
  70.    delete [nrows] vv;
  71. }
  72.  
  73. void CroutMatrix::lubksb(real* B, int mini)
  74. {
  75.    REPORT
  76.    int i,j; int ii=-1; real* ai=store;
  77.  
  78.    for (i=0;i<nrows;i++)
  79.    {
  80.       int ip=indx[i]; real sum=B[ip]; B[ip]=B[i];
  81.       if (ii>=0)
  82.       {
  83.          real* aij=ai+ii; real* bj=B+ii; j=i-ii;
  84.          while (j--) { sum -= *aij++ * *bj; bj++; }
  85.       }
  86.       else if (sum) ii=i;
  87.       B[i]=sum; ai += nrows;
  88.    }
  89.  
  90.    for (i=nrows-1;i>=mini;i--)
  91.    {
  92.       real* bj=B+i; ai -= nrows; real* ajx=ai+i; real sum=*bj; bj++;
  93.       j=nrows-i; while(--j) { sum -= *(++ajx) * *bj; bj++; }
  94.       B[i]=sum/(*(ai+i));
  95.    }
  96. }
  97.  
  98.  
  99. /****************************** scalar functions ****************************/
  100.  
  101. inline real square(real x) { return x*x; }
  102.  
  103. real GeneralMatrix::SumSquare()
  104. {
  105.    REPORT
  106.    real sum = 0.0; int i = storage; real* s = store;
  107.    while (i--) sum += square(*s++);
  108.    tDelete(); return sum;
  109. }
  110.  
  111. real GeneralMatrix::SumAbsoluteValue()
  112. {
  113.    REPORT
  114.    real sum = 0.0; int i = storage; real* s = store;
  115.    while (i--) sum += fabs(*s++);
  116.    tDelete(); return sum;
  117. }
  118.  
  119. real GeneralMatrix::MaximumAbsoluteValue()
  120. {
  121.    REPORT
  122.    real maxval = 0.0; int i = storage; real* s = store;
  123.    while (i--) { real a = fabs(*s++); if (maxval < a) maxval = a; }
  124.    tDelete(); return maxval;
  125. }
  126.  
  127. real SymmetricMatrix::SumSquare()
  128. {
  129.    REPORT
  130.    real sum1 = 0.0; real sum2 = 0.0; real* s = store; int nr = nrows;
  131.    for (int i = 0; i<nr; i++)
  132.    {
  133.       int j = i;
  134.       while (j--) sum2 += square(*s++);
  135.       sum1 += square(*s++);
  136.    }
  137.    tDelete(); return sum1 + 2.0 * sum2;
  138. }
  139.  
  140. real SymmetricMatrix::SumAbsoluteValue()
  141. {
  142.    REPORT
  143.    real sum1 = 0.0; real sum2 = 0.0; real* s = store; int nr = nrows;
  144.    for (int i = 0; i<nr; i++)
  145.    {
  146.       int j = i;
  147.       while (j--) sum2 += fabs(*s++);
  148.       sum1 += fabs(*s++);
  149.    }
  150.    tDelete(); return sum1 + 2.0 * sum2;
  151. }
  152.  
  153. real BaseMatrix::SumSquare()
  154. {
  155.    REPORT GeneralMatrix* gm = Evaluate();
  156.    real s = gm->SumSquare(); return s;
  157. }
  158.  
  159. real BaseMatrix::SumAbsoluteValue()
  160. {
  161.    REPORT GeneralMatrix* gm = Evaluate();
  162.    real s = gm->SumAbsoluteValue(); return s;
  163. }
  164.  
  165. real BaseMatrix::MaximumAbsoluteValue()
  166. {
  167.    REPORT GeneralMatrix* gm = Evaluate();
  168.    real s = gm->MaximumAbsoluteValue(); return s;
  169. }
  170.  
  171. real Matrix::Trace()
  172. {
  173.    REPORT
  174.    int i = nrows; int d = i+1;
  175.    if (i != ncols) MatrixError("trace of non-square matrix");
  176.    real sum = 0.0; real* s = store;
  177.    while (i--) { sum += *s; s += d; }
  178.    tDelete(); return sum;
  179. }
  180.  
  181. real DiagonalMatrix::Trace()
  182. {
  183.    REPORT
  184.    int i = nrows; real sum = 0.0; real* s = store;
  185.    while (i--) sum += *s++;
  186.    tDelete(); return sum;
  187. }
  188.  
  189. real SymmetricMatrix::Trace()
  190. {
  191.    REPORT
  192.    int i = nrows; real sum = 0.0; real* s = store; int j = 2;
  193.    while (i--) { sum += *s; s += j++; }
  194.    tDelete(); return sum;
  195. }
  196.  
  197. real LowerTriangularMatrix::Trace()
  198. {
  199.    REPORT
  200.    int i = nrows; real sum = 0.0; real* s = store; int j = 2;
  201.    while (i--) { sum += *s; s += j++; }
  202.    tDelete(); return sum;
  203. }
  204.  
  205. real UpperTriangularMatrix::Trace()
  206. {
  207.    REPORT
  208.    int i = nrows; real sum = 0.0; real* s = store;
  209.    while (i) { sum += *s; s += i--; }
  210.    tDelete(); return sum;
  211. }
  212.  
  213. real BaseMatrix::Trace()
  214. {
  215.    REPORT
  216.    GeneralMatrix* gm = Evaluate(MatrixType::Diag);
  217.    real sum = gm->Trace(); return sum;
  218. }
  219.  
  220. void LogAndSign::operator*=(real x)
  221. {
  222.    if (x > 0.0) { log_value += log(x); }
  223.    else if (x < 0.0) { log_value += log(-x); sign = -sign; }
  224.    else sign = 0;
  225. }
  226.  
  227. real LogAndSign::Value() { return sign * exp(log_value); }
  228.  
  229. LogAndSign DiagonalMatrix::LogDeterminant()
  230. {
  231.    REPORT
  232.    int i = nrows; LogAndSign sum; real* s = store;
  233.    while (i--) sum *= *s++;
  234.    tDelete(); return sum;
  235. }
  236.  
  237. LogAndSign LowerTriangularMatrix::LogDeterminant()
  238. {
  239.    REPORT
  240.    int i = nrows; LogAndSign sum; real* s = store; int j = 2;
  241.    while (i--) { sum *= *s; s += j++; }
  242.    tDelete(); return sum;
  243. }
  244.  
  245. LogAndSign UpperTriangularMatrix::LogDeterminant()
  246. {
  247.    REPORT
  248.    int i = nrows; LogAndSign sum; real* s = store;
  249.    while (i) { sum *= *s; s += i--; }
  250.    tDelete(); return sum;
  251. }
  252.  
  253. LogAndSign BaseMatrix::LogDeterminant()
  254. {
  255.    REPORT GeneralMatrix* gm = Evaluate();
  256.    LogAndSign sum = gm->LogDeterminant(); return sum;
  257. }
  258.  
  259. LogAndSign GeneralMatrix::LogDeterminant()
  260. {
  261.    REPORT
  262.    if (nrows != ncols) MatrixError("determinant of non-square matrix"); 
  263.    CroutMatrix C(*this); return C.LogDeterminant();
  264. }
  265.  
  266. LogAndSign CroutMatrix::LogDeterminant()
  267. {
  268.    REPORT
  269.    int i = nrows; int dd = i+1; LogAndSign sum; real* s = store;
  270.    while (i--) { sum *= *s; s += dd; }
  271.    if (!d) sum.ChangeSign(); return sum;
  272.  
  273. }
  274.