home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / newmat04 / part03 / newmat4.cxx < prev   
Encoding:
C/C++ Source or Header  |  1991-11-30  |  12.2 KB  |  404 lines
  1. //$$ newmat4.cxx       Constructors, ReDimension, basic utilities
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "newmat.hxx"
  8. #include "newmatrc.hxx"
  9.  
  10. //#define REPORT { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
  11.  
  12. #define REPORT {}
  13.  
  14. //#define REPORT1 { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
  15.  
  16. // REPORT1 constructors only - doesn't work in turbo and Borland C++
  17.  
  18. #define REPORT1 {}
  19.  
  20. //#define MONITOR(what,storage,store) \
  21. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  22.  
  23. #define MONITOR(what,store,storage) {}
  24.  
  25. /*************************** general utilities *************************/
  26.  
  27. static int tristore(int n)                      // els in triangular matrix
  28. { return (n*(n+1))/2; }
  29.  
  30.  
  31. /****************************** constructors ***************************/
  32.  
  33. GeneralMatrix::GeneralMatrix()
  34. { store=0; storage=0; nrows=0; ncols=0; tag=-1; }
  35.  
  36. GeneralMatrix::GeneralMatrix(int s)
  37. {
  38.    REPORT1
  39.    storage=s; tag=-1;
  40.    store = new real [storage]; MatrixErrorNoSpace(store);
  41.    MONITOR("Make (GenMatrix)",storage,store)
  42. }
  43.  
  44. Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
  45. { REPORT1 nrows=m; ncols=n; }
  46.  
  47. SymmetricMatrix::SymmetricMatrix(int n) : GeneralMatrix(tristore(n))
  48. { REPORT1 nrows=n; ncols=n; }
  49.  
  50. UpperTriangularMatrix::UpperTriangularMatrix(int n)
  51.    : GeneralMatrix(tristore(n))
  52. { REPORT1 nrows=n; ncols=n; }
  53.  
  54. LowerTriangularMatrix::LowerTriangularMatrix(int n)
  55.    : GeneralMatrix(tristore(n))
  56. { REPORT1 nrows=n; ncols=n; }
  57.  
  58. DiagonalMatrix::DiagonalMatrix(int m) : GeneralMatrix(m)
  59. { REPORT1 nrows=m; ncols=m; }
  60.  
  61. Matrix::Matrix(BaseMatrix& M)
  62. {
  63.    REPORT1 CheckConversion(M);
  64.    GeneralMatrix* gmx=M.Evaluate(MatrixType::Rect); GetMatrix(gmx);
  65. }
  66.  
  67. RowVector::RowVector(BaseMatrix& M) : Matrix(M)
  68. {
  69.    if (nrows!=1) MatrixError("Illegal conversion to row vector"); }
  70.  
  71. ColumnVector::ColumnVector(BaseMatrix& M) : Matrix(M)
  72. {
  73.    if (ncols!=1) MatrixError("Illegal conversion to column vector"); }
  74.  
  75. SymmetricMatrix::SymmetricMatrix(BaseMatrix& M)
  76. {
  77.    REPORT1 CheckConversion(M);
  78.    GeneralMatrix* gmx=M.Evaluate(MatrixType::Sym); GetMatrix(gmx);
  79. }
  80.  
  81. UpperTriangularMatrix::UpperTriangularMatrix(BaseMatrix& M)
  82. {
  83.    REPORT1 CheckConversion(M);
  84.    GeneralMatrix* gmx=M.Evaluate(MatrixType::UT); GetMatrix(gmx);
  85. }
  86.  
  87. LowerTriangularMatrix::LowerTriangularMatrix(BaseMatrix& M)
  88. {
  89.    REPORT1 CheckConversion(M);
  90.    GeneralMatrix* gmx=M.Evaluate(MatrixType::LT); GetMatrix(gmx);
  91. }
  92.  
  93. DiagonalMatrix::DiagonalMatrix(BaseMatrix& M)
  94. {
  95.    REPORT1 CheckConversion(M);
  96.    GeneralMatrix* gmx=M.Evaluate(MatrixType::Diag); GetMatrix(gmx);
  97. }
  98.  
  99. GeneralMatrix::~GeneralMatrix()
  100. {
  101.    if (store)
  102.    {
  103.       MONITOR("Free (GenMatrix)",storage,store)
  104.       delete [storage] store;
  105.    }
  106. }
  107.  
  108. CroutMatrix::CroutMatrix(BaseMatrix& m)
  109. {
  110.    REPORT1
  111.    GeneralMatrix* gm = m.Evaluate(MatrixType::Rect); GetMatrix(gm);
  112.    if (nrows!=ncols) MatrixError("Matrix not square");
  113.    d=TRUE; indx=new int [nrows]; MatrixErrorNoSpace(indx);
  114.    ludcmp();
  115. }
  116.  
  117. ReturnMatrix::ReturnMatrix(GeneralMatrix& gmx)
  118. {
  119.    REPORT1
  120.    gm = gmx.Type().New();  MatrixErrorNoSpace(gm);
  121.    gm->GetMatrix(&gmx);  gm->ReleaseAndDelete();
  122. }
  123.  
  124.  
  125. /**************************** ReDimension matrices ***************************/
  126.  
  127. void GeneralMatrix::ReDimension(int nr, int nc, int s)
  128. {
  129.    REPORT 
  130.    if (store)
  131.    {
  132.       MONITOR("Free (ReDimensi)",storage,store)
  133.       delete [storage] store;
  134.    }
  135.    storage=s; nrows=nr; ncols=nc; tag=-1;
  136.    store = new real [storage]; MatrixErrorNoSpace(store);
  137.    MONITOR("Make (ReDimensi)",storage,store)
  138. }
  139.  
  140. void Matrix::ReDimension(int nr, int nc)
  141. { REPORT GeneralMatrix::ReDimension(nr,nc,nr*nc); }
  142.  
  143. void SymmetricMatrix::ReDimension(int nr)
  144. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  145.  
  146. void UpperTriangularMatrix::ReDimension(int nr)
  147. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  148.  
  149. void LowerTriangularMatrix::ReDimension(int nr)
  150. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  151.  
  152. void DiagonalMatrix::ReDimension(int nr)
  153. { REPORT GeneralMatrix::ReDimension(nr,nr,nr); }
  154.  
  155. void RowVector::ReDimension(int nc)
  156. { REPORT GeneralMatrix::ReDimension(1,nc,nc); }
  157.  
  158. void ColumnVector::ReDimension(int nr)
  159. { REPORT GeneralMatrix::ReDimension(nr,1,nr); }
  160.  
  161.  
  162.  
  163. /********************* manipulate types, storage **************************/
  164.  
  165. int GeneralMatrix::search(const GeneralMatrix* s) const
  166. { REPORT return (s==this) ? 1 : 0; }
  167.  
  168. int AddedMatrix::search(const GeneralMatrix* s) const
  169. { REPORT return bm1->search(s) + bm2->search(s); }
  170.  
  171. int ShiftedMatrix::search(const GeneralMatrix* s) const
  172. { REPORT return bm->search(s); }
  173.  
  174. int NegatedMatrix::search(const GeneralMatrix* s) const
  175. { REPORT return bm->search(s); }
  176.  
  177. int ConstMatrix::search(const GeneralMatrix* s) const
  178. { REPORT return (s==cgm) ? 1 : 0; }
  179.  
  180. int ReturnMatrix::search(const GeneralMatrix* s) const
  181. { REPORT return (s==gm) ? 1 : 0; }
  182.  
  183. MatrixType Matrix::Type() const { return MatrixType::Rect; }
  184.  
  185. MatrixType SymmetricMatrix::Type() const { return MatrixType::Sym; }
  186.  
  187. MatrixType UpperTriangularMatrix::Type() const { return MatrixType::UT; }
  188.  
  189. MatrixType LowerTriangularMatrix::Type() const { return MatrixType::LT; }
  190.  
  191. MatrixType DiagonalMatrix::Type() const { return MatrixType::Diag; }
  192.  
  193. MatrixType RowVector::Type() const { return MatrixType::RowV; }
  194.  
  195. MatrixType ColumnVector::Type() const { return MatrixType::ColV; }
  196.  
  197. MatrixType CroutMatrix::Type() const { return MatrixType::Crout; }
  198.  
  199. MatrixType RowedMatrix::Type() const { return MatrixType::RowV; }
  200.  
  201. MatrixType ColedMatrix::Type() const { return MatrixType::ColV; }
  202.  
  203. MatrixType DiagedMatrix::Type() const { return MatrixType::Diag; }
  204.  
  205. MatrixType MatedMatrix::Type() const { return MatrixType::Rect; }
  206.  
  207. MatrixType GetSubMatrix::Type() const { return mt; }
  208.  
  209. MatrixType AddedMatrix::Type() const
  210. { REPORT return bm1->Type() + bm2->Type(); }
  211.  
  212. MatrixType MultipliedMatrix::Type() const
  213. { REPORT return bm1->Type() * bm2->Type(); }
  214.  
  215. MatrixType SolvedMatrix::Type() const
  216. { REPORT return bm1->Type().i() * bm2->Type(); }
  217.  
  218. MatrixType SubtractedMatrix::Type() const
  219. { REPORT return bm1->Type() - bm2->Type(); }
  220.  
  221. MatrixType ShiftedMatrix::Type() const
  222. { REPORT MatrixType mteqel(MatrixType::EqEl); return bm->Type() + mteqel; }
  223.  
  224. MatrixType ScaledMatrix::Type() const { REPORT return bm->Type(); }
  225. MatrixType TransposedMatrix::Type() const { REPORT return bm->Type().t(); }
  226. MatrixType NegatedMatrix::Type() const { REPORT return bm->Type(); }
  227. MatrixType InvertedMatrix::Type() const { REPORT return bm->Type().i(); }
  228. MatrixType ConstMatrix::Type() const { REPORT return cgm->Type(); }
  229. MatrixType ReturnMatrix::Type() const { REPORT return gm->Type(); }
  230.  
  231. int GeneralMatrix::NrowsV() const { return nrows; }
  232. int RowedMatrix::NrowsV() const { return 1; }
  233. int MatedMatrix::NrowsV() const { return nr; }
  234. int GetSubMatrix::NrowsV() const { return row_number; }
  235. int AddedMatrix::NrowsV() const  { return bm1->NrowsV(); }
  236. int ShiftedMatrix::NrowsV() const { return bm->NrowsV(); }
  237. int TransposedMatrix::NrowsV() const { return bm->NcolsV(); }
  238. int NegatedMatrix::NrowsV() const { return bm->NrowsV(); }
  239. int ColedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }
  240. int DiagedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }
  241. int ConstMatrix::NrowsV() const { return cgm->Nrows(); }
  242. int ReturnMatrix::NrowsV() const { return gm->Nrows(); }
  243.  
  244. int GeneralMatrix::NcolsV() const { return ncols; }
  245. int ColedMatrix::NcolsV() const { return 1; }
  246. int MatedMatrix::NcolsV() const { return nc; }
  247. int GetSubMatrix::NcolsV() const { return col_number; }
  248. int AddedMatrix::NcolsV() const  { return bm2->NcolsV(); }
  249. int ShiftedMatrix::NcolsV() const { return bm->NcolsV(); }
  250. int TransposedMatrix::NcolsV() const { return bm->NrowsV(); }
  251. int NegatedMatrix::NcolsV() const { return bm->NcolsV(); }
  252. int RowedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
  253. int DiagedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
  254. int ConstMatrix::NcolsV() const { return cgm->Ncols(); }
  255. int ReturnMatrix::NcolsV() const { return gm->Ncols(); }
  256.  
  257.  
  258. //  Rules regarding tDelete, reuse, GetStore
  259. //    All matrices processed during expression evaluation must be subject
  260. //    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
  261. //    If reuse returns TRUE the matrix must be reused.
  262. //    GetMatrix(gm) always calls gm->GetStore()
  263. //    gm->Evaluate obeys rules
  264. //    bm->Evaluate obeys rules for matrices in bm structure
  265.  
  266. void GeneralMatrix::tDelete()
  267. {
  268.    if (tag<0)
  269.    {
  270.       if (tag<-1) { REPORT store=0; delete this; return; }  // borrowed
  271.       else { REPORT return; }
  272.    }
  273.    if (tag==1)
  274.    {
  275.       REPORT  MONITOR("Free   (tDelete)",storage,store)
  276.       if (store) delete [storage] store; store=0; tag=-1; return;
  277.    }
  278.    if (tag==0) { REPORT delete this; return; }
  279.    REPORT tag--; return;
  280. }
  281.  
  282. BOOL GeneralMatrix::reuse()
  283. {
  284.    if (tag<-1)
  285.    {
  286.       REPORT
  287.       real* s = new real [storage]; MatrixErrorNoSpace(s);
  288.       MONITOR("Make     (reuse)",storage,s)
  289.       real* s1=store; real* s2=s; int i=storage;
  290.       while(i--) *s2++ = *s1++;
  291.       store=s; tag=0; return TRUE;
  292.    }
  293.    if (tag<0) { REPORT return FALSE; }
  294.    if (tag<=1)  { REPORT return TRUE; }
  295.    REPORT tag--; return FALSE;
  296. }
  297.  
  298. real* GeneralMatrix::GetStore()
  299. {
  300.    if (tag<0 || tag>1)
  301.    {
  302.       real* s = new real [storage]; MatrixErrorNoSpace(s);
  303.       MONITOR("Make  (GetStore)",storage,s)
  304.       real* s1=store; real* s2=s; int i=storage;
  305.       while(i--) *s2++ = *s1++;
  306.       if (tag>1) { REPORT tag--; }
  307.       else if (tag < -1) { REPORT store=0; delete this; } // borrowed store
  308.       else { REPORT }
  309.       return s;
  310.    }
  311.    real* s=store; store=0;
  312.    if (tag==0) { REPORT delete this; }
  313.    else { REPORT tag=-1; }
  314.    return s;
  315. }
  316.  
  317. #ifndef __ZTC__
  318. void GeneralMatrix::GetMatrixC(const GeneralMatrix* gmx)
  319. {
  320.    REPORT tag=-1;
  321.    nrows=gmx->nrows; ncols=gmx->ncols; storage=gmx->storage;
  322.  
  323.    store = new real [storage]; MatrixErrorNoSpace(store);
  324.    MONITOR("Make (GetMatrix)",storage,store)
  325.    real* s1=gmx->store; real* s2=store; int i=storage;
  326.    while(i--) *s2++ = *s1++;
  327. }
  328. #endif
  329.  
  330. void GeneralMatrix::GetMatrix(GeneralMatrix* gmx)
  331. {
  332.    REPORT  tag=-1; nrows=gmx->Nrows(); ncols=gmx->Ncols();
  333.    storage=gmx->storage; store=gmx->GetStore(); 
  334. }
  335.  
  336. GeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
  337. // Copy storage of *this to storage of *gmx. Then convert to type mt.
  338. // If mt == 0 just let *gm point to storage of *this if tag<0.
  339. {
  340.    if ((int)mt==0)
  341.    {
  342.       if (tag == -1) { REPORT gmx->tag = -2; gmx->store = store; }
  343.       else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
  344.    }
  345.    else if (mt!=gmx->Type())
  346.    {
  347.       REPORT gmx->tag = -2; gmx->store = store;
  348.       gmx = gmx->Evaluate(mt); gmx->tag = 0; tDelete();
  349.    }
  350.    else { REPORT  gmx->tag = 0; gmx->store = GetStore(); }
  351.  
  352.    return gmx;
  353. }   
  354.  
  355. void GeneralMatrix::Eq(BaseMatrix& X, MatrixType mt)
  356. // Count number of references to this in X.
  357. // If zero delete storage in X;
  358. // otherwise tag X to show when storage can be deleted
  359. // evaluate X and copy to current object
  360. {
  361.    int counter=X.search(this);
  362.    if (counter==0)
  363.    {
  364.       REPORT  MONITOR("Free (operator=)",storage,store)
  365.       if (store) { REPORT delete [storage] store; storage=0; }
  366.    }
  367.    else { REPORT Release(counter); }
  368.    GeneralMatrix* gmx = X.Evaluate(mt);
  369.    if (gmx!=this) { REPORT GetMatrix(gmx); }
  370.    else { REPORT }
  371.    Protect();
  372. }
  373.  
  374. void GeneralMatrix::Inject(const GeneralMatrix& X)
  375. // copy stored values of X; otherwise leave els of *this unchanged
  376. {
  377.    REPORT
  378.    CheckConversion(X);
  379.    if (nrows != X.nrows || ncols != X.ncols)
  380.       MatrixError("Inject: dimensions don't match");
  381.    MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
  382.    MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
  383.    int i=nrows;
  384.    while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
  385. }  
  386.  
  387. void GeneralMatrix::CheckConversion(const BaseMatrix& M)
  388. { if (!(this->Type() >= M.Type())) MatrixError("Illegal Conversion"); }
  389.  
  390. /************************* nricMatrix routines *****************************/
  391.  
  392. void nricMatrix::MakeRowPointer()
  393. {
  394.    row_pointer = new real* [nrows]; MatrixErrorNoSpace(row_pointer);
  395.    real* s = Store() - 1; int i = nrows; real** rp = row_pointer;
  396.    while (i--) { *rp++ = s; s+=ncols; }
  397. }
  398.  
  399. void nricMatrix::DeleteRowPointer()
  400. { if (nrows) delete [nrows] row_pointer; }
  401.  
  402. void GeneralMatrix::CheckStore() const
  403. { if (!store) MatrixError("NRIC accessing matrix with dimensions not set"); }
  404.