home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / newmat04 / part03 / newmat1.cxx next >
Encoding:
C/C++ Source or Header  |  1991-11-30  |  7.6 KB  |  283 lines
  1. //$$ newmat1.cxx   Utilities and matrix type functions
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #define WANT_STREAM                     // to include stream package
  6.  
  7. #include "include.hxx"
  8.  
  9. #include "newmat.hxx"
  10.  
  11.  
  12. /**************************** error handlers *******************************/
  13.  
  14. void MatrixError(char* message)
  15. {
  16.    cerr << "\nError detected by matrix package: " << message << "\n";
  17.    exit(0);
  18. }
  19.  
  20. void MatrixErrorNoSpace(void* v)
  21. // give error message if v is null
  22. {
  23.    if (v) return;
  24.    cerr << "\nError detected by matrix package: no space\n";
  25.    exit(0);
  26. }
  27.  
  28. /************************* ExeCounter functions *****************************/
  29.  
  30.  
  31. // the next statment may need to be deleted form some compilers
  32.  
  33. int ExeCounter::nreports = 0;
  34.  
  35. ExeCounter::ExeCounter(int xl, int xf) : line(xl), fileid(xf), nexe(0) {}
  36.  
  37. ExeCounter::~ExeCounter()
  38. {
  39.    nreports++;
  40.    cout << nreports << "  " << fileid << "  " << line << "  " << nexe << "\n";
  41. }
  42.  
  43.  
  44. /************************* MatrixType functions *****************************/
  45.  
  46.  
  47. BOOL MatrixType::operator==(const MatrixType& t) const
  48. { return (type == t.type); }
  49.  
  50. BOOL MatrixType::operator!=(const MatrixType& t) const
  51. { return (type != t.type); }
  52.  
  53. MatrixType MatrixType::operator+(const MatrixType& mt) const
  54. {
  55.    Type type1=mt.type;
  56.    if (type==UnSp || type1==UnSp) return UnSp;
  57.    if (type==type1) return type;   // won't work for orthogonal
  58.    if (type==ColV || type1==ColV) return ColV;
  59.    if (type==RowV || type1==RowV) return RowV;
  60.    if (type==EqEl && (type1==Sym || type1==Diag)) return Sym;
  61.    if (type1==EqEl && (type==Sym || type==Diag)) return Sym;
  62.    if (type==EqEl || type1==EqEl) return Rect;
  63.    if (type==Diag) return type1;   // won't work for assym
  64.    if (type1==Diag) return type;
  65.    return Rect;
  66. }
  67.  
  68. MatrixType MatrixType::operator-(const MatrixType& mt) const
  69. {
  70.    return *this+mt;        // won't work for orthogonal or pos def
  71. }
  72.  
  73. MatrixType MatrixType::operator*(const MatrixType& mt) const
  74. {
  75.    Type type1=mt.type;
  76.    if (type==UnSp || type1==UnSp) return UnSp;
  77.    if (type1==ColV) { if (type==RowV) return Diag; else return ColV; }
  78.    if (type==RowV) return RowV;
  79.    if (type==Sym || type1==Sym) return Rect;
  80.    if (type==type1) return type;    // won't work for pos def or assym
  81.    if (type==EqEl || type1==EqEl) return Rect;
  82.    if (type==Diag) return type1;
  83.    if (type1==Diag) return type;
  84.    return Rect;
  85. }
  86.  
  87. BOOL MatrixType::operator>=(const MatrixType& mt) const
  88. {
  89.    Type type1=mt.type;
  90.    if (type==type1) return TRUE;    // by definition
  91.    if (type==Rect || type==RowV || type==ColV) return TRUE;
  92.    if (type==EqEl) return FALSE;
  93.    if (type1==Diag) return TRUE;    // won't work for pos def or assym
  94.    return FALSE;
  95. }
  96.  
  97. MatrixType MatrixType::i() const
  98. {
  99.    switch (type)
  100.    {
  101.    case UnSp:  return UnSp;
  102.    case UT:    return UT;
  103.    case LT:    return LT;
  104.    case Rect:  return Rect;
  105.    case Sym:   return Sym;
  106.    case Diag:  return Diag;
  107.    case RowV:  return Diag;
  108.    case ColV:  return Diag;
  109.    case EqEl:  return Diag;
  110.    case Crout: return UnSp;
  111.    }
  112. }
  113.  
  114. MatrixType MatrixType::t() const
  115. {
  116.    switch (type)
  117.    {
  118.    case UnSp:  return UnSp;
  119.    case UT:    return LT;
  120.    case LT:    return UT;
  121.    case Rect:  return Rect;
  122.    case Sym:   return Sym;
  123.    case Diag:  return Diag;
  124.    case RowV:  return ColV;
  125.    case ColV:  return RowV;
  126.    case EqEl:  return EqEl;
  127.    case Crout: return UnSp;
  128.    }
  129. }
  130.  
  131. MatrixType MatrixType::sub() const
  132. {
  133.    switch (type)
  134.    {
  135.    case UnSp:  return UnSp;
  136.    case UT:    return Rect;
  137.    case LT:    return Rect;
  138.    case Rect:  return Rect;
  139.    case Sym:   return Rect;
  140.    case Diag:  return Rect;
  141.    case RowV:  return RowV;
  142.    case ColV:  return ColV;
  143.    case EqEl:  return EqEl;
  144.    case Crout: return UnSp;
  145.    }
  146. }
  147.  
  148. MatrixType MatrixType::ssub() const
  149. {
  150.    return *this;                  // won't work for selection matrix
  151. }
  152.  
  153. MatrixType::operator char*() const
  154. {
  155.    switch (type)
  156.    {
  157.    case UnSp:  return "UnSp ";
  158.    case UT:    return "UT   ";
  159.    case LT:    return "LT   ";
  160.    case Rect:  return "Rect ";
  161.    case Sym:   return "Sym  ";
  162.    case Diag:  return "Diag ";
  163.    case RowV:  return "RowV ";
  164.    case ColV:  return "ColV ";
  165.    case EqEl:  return "EqEl ";
  166.    case Crout: return "Crout";
  167.    }
  168. }
  169.  
  170. GeneralMatrix* MatrixType::New() const
  171. {
  172.    GeneralMatrix* gm;
  173.    switch (type)
  174.    {
  175.    case UnSp:  MatrixError("Can't build matrix type UnSp");
  176.    case UT:    gm = new UpperTriangularMatrix; break;
  177.    case LT:    gm = new LowerTriangularMatrix; break;
  178.    case Rect:  gm = new Matrix; break;
  179.    case Sym:   gm = new SymmetricMatrix; break;
  180.    case Diag:  gm = new DiagonalMatrix; break;
  181.    case RowV:  gm = new RowVector; break;
  182.    case ColV:  gm = new ColumnVector; break;
  183.    case EqEl:  MatrixError("Can't build matrix type EqEl");
  184.    case Crout: MatrixError("Can't build matrix type Crout");
  185.    }
  186.    MatrixErrorNoSpace(gm);
  187.    gm->Protect(); return gm;
  188. }
  189.  
  190. GeneralMatrix* MatrixType::New(int nr, int nc) const
  191. {
  192.    GeneralMatrix* gm;
  193.    switch (type)
  194.    {
  195.    case UnSp:  MatrixError("Can't build matrix type UnSp");
  196.    case UT:    gm = new UpperTriangularMatrix(nr); break;
  197.    case LT:    gm = new LowerTriangularMatrix(nr); break;
  198.    case Rect:  gm = new Matrix(nr, nc); break;
  199.    case Sym:   gm = new SymmetricMatrix(nr); break;
  200.    case Diag:  gm = new DiagonalMatrix(nr); break;
  201.    case RowV:  if (nr!=1) MatrixError("Illegal Row Vector");
  202.            gm = new RowVector(nc); break;
  203.    case ColV:  if (nc!=1) MatrixError("Illegal Column Vector");
  204.            gm = new ColumnVector(nr); break;
  205.    case EqEl:  MatrixError("Can't build matrix type EqEl");
  206.    case Crout: MatrixError("Can't build matrix type Crout");
  207.    }
  208.    MatrixErrorNoSpace(gm);
  209.    gm->Protect(); return gm;
  210. }
  211.  
  212. void TestTypeAdd()
  213. {
  214.    MatrixType list[] = { MatrixType::UT,
  215.                          MatrixType::LT,
  216.                          MatrixType::Rect,
  217.                          MatrixType::Sym,
  218.                          MatrixType::Diag,
  219.                          MatrixType::RowV,
  220.                          MatrixType::ColV,
  221.              MatrixType::EqEl };
  222.  
  223.    cout << "+     ";
  224.    for (int i=0; i<MatrixType::nTypes(); i++) cout << (char*)list[i] << " ";
  225.    cout << "\n";
  226.    for (i=0; i<MatrixType::nTypes(); i++)
  227.    {
  228.       cout << (char*)(list[i]) << " ";
  229.       for (int j=0; j<MatrixType::nTypes(); j++)
  230.      cout << (char*)(list[j]+list[i]) << " ";
  231.       cout << "\n";
  232.    }
  233.    cout << "\n";
  234. }
  235.  
  236. void TestTypeMult()
  237. {
  238.    MatrixType list[] = { MatrixType::UT,
  239.                          MatrixType::LT,
  240.                          MatrixType::Rect,
  241.                          MatrixType::Sym,
  242.                          MatrixType::Diag,
  243.                          MatrixType::RowV,
  244.                          MatrixType::ColV,
  245.              MatrixType::EqEl };
  246.    cout << "*     ";
  247.    for (int i=0; i<MatrixType::nTypes(); i++)
  248.       cout << (char*)list[i] << " ";
  249.    cout << "\n";
  250.    for (i=0; i<MatrixType::nTypes(); i++)
  251.    {
  252.       cout << (char*)list[i] << " ";
  253.       for (int j=0; j<MatrixType::nTypes(); j++)
  254.      cout << (char*)(list[j]*list[i]) << " ";
  255.       cout << "\n";
  256.    }
  257.    cout << "\n";
  258. }
  259.  
  260. void TestTypeOrder()
  261. {
  262.    MatrixType list[] = { MatrixType::UT,
  263.                          MatrixType::LT,
  264.                          MatrixType::Rect,
  265.                          MatrixType::Sym,
  266.                          MatrixType::Diag,
  267.                          MatrixType::RowV,
  268.                          MatrixType::ColV,
  269.              MatrixType::EqEl };
  270.    cout << ">=    ";
  271.    for (int i = 0; i<MatrixType::nTypes(); i++)
  272.       cout << (char*)list[i] << " ";
  273.    cout << "\n";
  274.    for (i=0; i<MatrixType::nTypes(); i++)
  275.    {
  276.       cout << (char*)list[i] << " ";
  277.       for (int j=0; j<MatrixType::nTypes(); j++)
  278.      cout << ((list[j]>=list[i]) ? "Yes   " : "No    ");
  279.       cout << "\n";
  280.    }
  281.    cout << "\n";
  282. }
  283.