home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / BANDMAT.CXX next >
C/C++ Source or Header  |  1992-11-29  |  10KB  |  360 lines
  1. //$$ bandmat.cxx                     Band matrix definitions
  2.  
  3. // Copyright (C) 1991,2: R B Davies
  4.  
  5. #define WANT_MATH                    // include.h will get math fns
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmat.h"
  10. #include "newmatrc.h"
  11.  
  12. //#define REPORT { static ExeCounter ExeCount(__LINE__,4); ++ExeCount; }
  13.  
  14. #define REPORT {}
  15.  
  16. //#define REPORT1 { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
  17.  
  18. // REPORT1 constructors only - doesn't work in turbo and Borland C++
  19.  
  20. #define REPORT1 {}
  21.  
  22. //#define MONITOR(what,storage,store) \
  23. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  24.  
  25. #define MONITOR(what,store,storage) {}
  26.  
  27. BandMatrix::BandMatrix(const BaseMatrix& M)
  28. {
  29.    REPORT1 CheckConversion(M);
  30.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::BM);
  31.    GetMatrix(gmx); CornerClear();
  32. }
  33.  
  34. void BandMatrix::SetParameters(const GeneralMatrix* gmx)
  35. {
  36.    MatrixBandWidth bw = gmx->BandWidth();
  37.    lower = bw.lower; upper = bw.upper;
  38. }
  39.  
  40. void BandMatrix::ReDimension(int n, int lb, int ub)
  41. {
  42.    REPORT
  43.    Tracer tr("BandMatrix::ReDimension");
  44.    if (lb<0 || ub<0) Throw(ProgramException("Undefined bandwidth"));
  45.    lower = (lb<=n) ? lb : n-1; upper = (ub<=n) ? ub : n-1;
  46.    GeneralMatrix::ReDimension(n,n,n*(lower+1+upper)); CornerClear();
  47. }
  48.  
  49. void BandMatrix::operator=(const BaseMatrix& X)
  50. { REPORT CheckConversion(X); Eq(X,MatrixType::BM); CornerClear(); }
  51.  
  52. void BandMatrix::CornerClear() const
  53. {
  54.    // set unused parts of BandMatrix to zero
  55.    REPORT
  56.    int i = lower; Real* s = store; int bw = lower + 1 + upper;
  57.    while (i)
  58.       { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
  59.    i = upper; s = store + storage;
  60.    while (i)
  61.       { int j = i--; Real* sj = s; s -= bw; while (j--) *(--sj) = 0.0; }
  62. }
  63.  
  64. MatrixBandWidth MatrixBandWidth::operator+(const MatrixBandWidth& bw) const
  65. {
  66.    int l = bw.lower; int u = bw.upper;
  67.    l = (lower < 0 || l < 0) ? -1 : (lower > l) ? lower : l;
  68.    u = (upper < 0 || u < 0) ? -1 : (upper > u) ? upper : u;
  69.    return MatrixBandWidth(l,u);
  70. }
  71.  
  72. MatrixBandWidth MatrixBandWidth::operator*(const MatrixBandWidth& bw) const
  73. {
  74.    int l = bw.lower; int u = bw.upper;
  75.    l = (lower < 0 || l < 0) ? -1 : lower+l;
  76.    u = (upper < 0 || u < 0) ? -1 : upper+u;
  77.    return MatrixBandWidth(l,u);
  78. }
  79.  
  80. UpperBandMatrix::UpperBandMatrix(const BaseMatrix& M)
  81. {
  82.    REPORT1 CheckConversion(M);
  83.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UB);
  84.    GetMatrix(gmx); CornerClear();
  85. }
  86.  
  87. void UpperBandMatrix::operator=(const BaseMatrix& X)
  88. { REPORT CheckConversion(X); Eq(X,MatrixType::UB); CornerClear(); }
  89.  
  90. LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M)
  91. {
  92.    REPORT1 CheckConversion(M);
  93.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LB);
  94.    GetMatrix(gmx); CornerClear();
  95. }
  96.  
  97. void LowerBandMatrix::operator=(const BaseMatrix& X)
  98. { REPORT CheckConversion(X); Eq(X,MatrixType::LB); CornerClear(); }
  99.  
  100. BandLUMatrix::BandLUMatrix(const BaseMatrix& m)
  101. {
  102.    REPORT1
  103.    Tracer tr("BandLUMatrix");
  104.    GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::BM);
  105.    GetMatrix(gm);
  106.    m1 = ((BandMatrix*)gm)->lower; m2 = ((BandMatrix*)gm)->upper;
  107.    if (nrows!=ncols) Throw(NotSquareException(*this));
  108.    d = TRUE; sing = FALSE;
  109.    indx = new int [nrows]; MatrixErrorNoSpace(indx);
  110.    MONITOR_INT_NEW("Index (BndLUMat)",nrows,indx)
  111.    storage2 = nrows * m1;
  112.    store2 = new Real [storage2]; MatrixErrorNoSpace(store2);
  113.    MONITOR_REAL_NEW("Make (BandLUMat)",storage2,store2)
  114.    ludcmp();
  115. }
  116.  
  117. BandLUMatrix::~BandLUMatrix()
  118. {
  119.    MONITOR_INT_DELETE("Index (BndLUMat)",nrows,indx)
  120.    MONITOR_REAL_DELETE("Delete (BndLUMt)",storage2,store2)
  121. #ifdef Version21
  122.    delete [] indx; delete [] store2;
  123. #else
  124.    delete [nrows] indx; delete [storage2] store2;
  125. #endif
  126. }
  127.  
  128. MatrixType BandLUMatrix::Type() const { return MatrixType::BC; }
  129.  
  130.  
  131. LogAndSign BandLUMatrix::LogDeterminant() const
  132. {
  133.    if (sing) return 0.0;
  134.    Real* a = store; int w = m1+1+m2; LogAndSign sum; int i = nrows;
  135.    while (i--) { sum *= *a; a += w; }
  136.    if (!d) sum.ChangeSign(); return sum;
  137. }
  138.  
  139. GeneralMatrix* BandMatrix::MakeSolver()
  140. {
  141.    REPORT
  142.    GeneralMatrix* gm = new BandLUMatrix(*this);
  143.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  144. }
  145.  
  146.  
  147. void BandLUMatrix::ludcmp()
  148. {
  149.    REPORT
  150.    Real* a = store;
  151.    int i = m1; int j = m2; int k; int n = nrows; int w = m1 + 1 + m2;
  152.    while (i)
  153.    {
  154.       Real* ai = a + i;
  155.       k = ++j; while (k--) *a++ = *ai++;
  156.       k = i--; while (k--) *a++ = 0.0;
  157.    }
  158.  
  159.    a = store; int l = m1;
  160.    for (k=0; k<n; k++)
  161.    {
  162.       Real x = *a; i = k; Real* aj = a;
  163.       if (l < n) l++;
  164.       for (j=k+1; j<l; j++)
  165.          { aj += w; if (fabs(x) < fabs(*aj)) { x = *aj; i = j; } }
  166.       indx[k] = i;
  167.       if (x==0) { sing = TRUE; return; }
  168.       if (i!=k)
  169.       {
  170.          d = !d; Real* ak = a; Real* ai = store + i * w; j = w;
  171.          while (j--) { x = *ak; *ak++ = *ai; *ai++ = x; }
  172.       }
  173.       aj = a + w; Real* m = store2 + m1 * k;
  174.       for (j=k+1; j<l; j++)
  175.       {
  176.          *m++ = x = *aj / *a; i = w; Real* ak = a;
  177.      while (--i) { Real* aj1 = aj++; *aj1 = *aj - x * *(++ak); }
  178.          *aj++ = 0.0;
  179.       }
  180.       a += w;
  181.    }
  182. }
  183.  
  184. void BandLUMatrix::lubksb(Real* B, int mini)
  185. {
  186.    REPORT
  187.    Tracer tr("BandLUMatrix::lubksb");
  188.    if (sing) Throw(SingularException(*this));
  189.    int n = nrows; int l = m1; int w = m1 + 1 + m2;
  190.  
  191.    for (int k=0; k<n; k++)
  192.    {
  193.       int i = indx[k];
  194.       if (i!=k) { Real x=B[k]; B[k]=B[i]; B[i]=x; }
  195.       if (l<n) l++;
  196.       Real* m = store2 + k*m1; Real* b = B+k; Real* bi = b;
  197.       for (i=k+1; i<l; i++)  *(++bi) -= *m++ * *b;
  198.    }
  199.  
  200.    l = -m1;
  201.    for (int i = n-1; i>=mini; i--)
  202.    {
  203.       Real* b = B + i; Real* bk = b; Real x = *bk;
  204.       Real* a = store + w*i; Real y = *a;
  205.       int k = l+m1; while (k--) x -=  *(++a) * *(++bk);
  206.       *b = x / y;
  207.       if (l < m2) l++;
  208.    }
  209. }
  210.  
  211. void BandLUMatrix::Solver(MatrixRowCol& mcout, const MatrixRowCol& mcin)
  212. {
  213.    REPORT
  214.    Real* el = mcin.store; int i = mcin.skip;
  215.    while (i--) *el++ = 0.0;
  216.    el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
  217.    while (i--) *el++ = 0.0;
  218.    lubksb(mcin.store, mcout.skip);
  219. }
  220.  
  221. // Do we need check for entirely zero output?
  222.  
  223.  
  224. void UpperBandMatrix::Solver(MatrixRowCol& mcout,
  225.    const MatrixRowCol& mcin)
  226. {
  227.    REPORT
  228.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  229.    while (i-- > 0) *elx++ = 0.0;
  230.    int nr = mcin.skip+mcin.storage; elx = mcin.store+nr; Real* el = elx;
  231.    int j = mcout.skip+mcout.storage-nr; i = nr-mcout.skip;
  232.    while (j-- > 0) *elx++ = 0.0;
  233.  
  234.    Real* Ael = store + (upper+1)*(i-1)+1; j = 0;
  235.    while (i-- > 0)
  236.    {
  237.       elx = el; Real sum = 0.0; int jx = j;
  238.       while (jx--) sum += *(--Ael) * *(--elx);
  239.       elx--; *elx = (*elx - sum) / *(--Ael);
  240.       if (j<upper) Ael -= upper - (++j); else el--;
  241.    }
  242. }
  243.  
  244. void LowerBandMatrix::Solver(MatrixRowCol& mcout,
  245.    const MatrixRowCol& mcin)
  246. {
  247.    REPORT
  248.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  249.    while (i-- > 0) *elx++ = 0.0;
  250.    int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.store+i;
  251.    int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
  252.    while (j-- > 0) *elx++ = 0.0;
  253.  
  254.    Real* el = mcin.store+nc; Real* Ael = store + (lower+1)*nc + lower; j = 0;
  255.    while (i-- > 0)
  256.    {
  257.       elx = el; Real sum = 0.0; int jx = j;
  258.       while (jx--) sum += *Ael++ * *elx++;
  259.       *elx = (*elx - sum) / *Ael++;
  260.       if (j<lower) Ael += lower - (++j); else el++;
  261.    }
  262. }
  263.  
  264.  
  265. LogAndSign BandMatrix::LogDeterminant() const
  266. {
  267.    REPORT
  268.    BandLUMatrix C(*this); return C.LogDeterminant();
  269. }
  270.  
  271. LogAndSign LowerBandMatrix::LogDeterminant() const
  272. {
  273.    REPORT
  274.    int i = nrows; LogAndSign sum; Real* s = store + lower; int j = lower + 1;
  275.    while (i--) { sum *= *s; s += j; }
  276.    ((GeneralMatrix&)*this).tDelete(); return sum;
  277. }
  278.  
  279. LogAndSign UpperBandMatrix::LogDeterminant() const
  280. {
  281.    REPORT
  282.    int i = nrows; LogAndSign sum; Real* s = store; int j = upper + 1;
  283.    while (i--) { sum *= *s; s += j; }
  284.    ((GeneralMatrix&)*this).tDelete(); return sum;
  285. }
  286.  
  287. GeneralMatrix* SymmetricBandMatrix::MakeSolver()
  288. {
  289.    REPORT
  290.    GeneralMatrix* gm = new BandLUMatrix(*this);
  291.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  292. }
  293.  
  294. SymmetricBandMatrix::SymmetricBandMatrix(const BaseMatrix& M)
  295. {
  296.    REPORT1 CheckConversion(M);
  297.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::SB);
  298.    GetMatrix(gmx);
  299. }
  300.  
  301. GeneralMatrix* SymmetricBandMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  302. { REPORT  return Evaluate(mt); }
  303.  
  304. LogAndSign SymmetricBandMatrix::LogDeterminant() const
  305. {
  306.    REPORT
  307.    BandLUMatrix C(*this); return C.LogDeterminant();
  308. }
  309.  
  310. void SymmetricBandMatrix::SetParameters(const GeneralMatrix* gmx)
  311. { lower = gmx->BandWidth().lower; }
  312.  
  313. void SymmetricBandMatrix::ReDimension(int n, int lb)
  314. {
  315.    REPORT
  316.    Tracer tr("SymmetricBandMatrix::ReDimension");
  317.    if (lb<0) Throw(ProgramException("Undefined bandwidth"));
  318.    lower = (lb<=n) ? lb : n-1;
  319.    GeneralMatrix::ReDimension(n,n,n*(lower+1));
  320. }
  321.  
  322. void SymmetricBandMatrix::operator=(const BaseMatrix& X)
  323. { REPORT CheckConversion(X); Eq(X,MatrixType::SB); }
  324.  
  325. void SymmetricBandMatrix::CornerClear() const
  326. {
  327.    // set unused parts of BandMatrix to zero
  328.    REPORT
  329.    int i = lower; Real* s = store; int bw = lower + 1;
  330.    while (i)
  331.       { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
  332. }
  333.  
  334. MatrixBandWidth SymmetricBandMatrix::BandWidth() const
  335.    { return MatrixBandWidth(lower,lower); }
  336.  
  337. inline Real square(Real x) { return x*x; }
  338.  
  339.  
  340. Real SymmetricBandMatrix::SumSquare() const
  341. {
  342.    REPORT
  343.    CornerClear();
  344.    Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
  345.    while (i--)
  346.       { int j = l; while (j--) sum2 += square(*s++); sum1 += square(*s++); }
  347.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  348. }
  349.  
  350. Real SymmetricBandMatrix::SumAbsoluteValue() const
  351. {
  352.    REPORT
  353.    CornerClear();
  354.    Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
  355.    while (i--)
  356.       { int j = l; while (j--) sum2 += fabs(*s++); sum1 += fabs(*s++); }
  357.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  358. }
  359.  
  360.