home *** CD-ROM | disk | FTP | other *** search
/ Power Programming / powerprogramming1994.iso / progtool / microcrn / issue_40.arc / DAIMS.ARC / MATRIX1.CPP < prev    next >
Text File  |  1988-02-10  |  15KB  |  630 lines
  1. #include "matrix.hxx"
  2. #include <math.h>
  3. #include <ctype.h>
  4. #include <time.h>
  5. #ifndef DOS
  6. #include <sys/param.h>
  7. #endif
  8. #define UPCASE(X)  (X = islower(X) ? toupper(X) : X) /* force a char to upper case */
  9.  
  10. /* 
  11. -*++ matrix::matrix(): contructs a matrix of any size and inits to 0
  12. ** 
  13. ** (*++ history: 
  14. **      4 Dec 87    Bruce Eckel    Creation date
  15. ** ++*)
  16. ** 
  17. ** (*++ detailed: 
  18. ** ++*)
  19. */
  20.  
  21.   matrix::matrix(int rows = 1 , int cols = 1, double initval = 0)
  22. { 
  23.   p = new matrep;
  24.   p->m_vec = new _Vector*[rows];
  25.   p->m_rows = rows;
  26.   p->m_cols = cols;
  27.   for (int i = 0; i < rows; i++)
  28.     p->m_vec[i] = new _Vector(cols);
  29.   for (i=0; i < rows; i++)
  30.     for (int j = 0; j < cols; j++)
  31.       (*this)[i][j] = initval;
  32.   p->n = 1;            /* so far, only one reference to this data */
  33. }
  34.  
  35.  
  36. /* 
  37. -*++ matrix::matrix(): adds a new reference to an existing matrix
  38. ** 
  39. ** (*++ history: 
  40. **      4 Dec 87    Bruce Eckel    Creation date
  41. ** ++*)
  42. ** 
  43. ** (*++ detailed: "Reference counting" is a method of managing garbage.
  44. ** Instead of creating a whole new object every time you use (for instance) an
  45. ** equal sign, it makes more sense just to point at the old object.  This
  46. ** means you can have several objects all pointing at the same thing.  When
  47. ** you start destroying objects, you don't want to destroy the original thing
  48. ** until you only have one reference left.
  49. ** ++*)
  50. */
  51.  
  52.   matrix::matrix(matrix& x)
  53. {
  54.   x.p->n++;    /* we're adding another reference, so increase */
  55.   p = x.p;    /* the count. */
  56. }
  57.  
  58.  
  59. /* 
  60. -*++ ch_find(): search an input stream for a char, ignoring comments
  61. ** 
  62. ** (*++ history: 
  63. **      4 Dec 87    Bruce Eckel    Creation date
  64. ** ++*)
  65. ** 
  66. ** (*++ detailed: searches an input stream for c, ignoring comments starting
  67. ** with '#' to end of line.  Errfunc is the address of an error function,
  68. ** fname is  the name of the source file. 
  69. ** ++*)
  70. */
  71.  
  72.     void ch_find(istream& from, char ch, void (*errfunc)(char *), char * fname)
  73. {
  74.   char c;
  75.   UPCASE(ch);
  76.   while(TRUE){        
  77.     while(from >> c, UPCASE(c), c != ch && c != '#') 
  78.       if(from.eof()) errfunc(fname); /* find ch or comment start */
  79.     if (c == ch) break;
  80.     /* If we made it here, a '#' comment start was found */
  81.     while (from.get(c), c != EOL) /* 'get' doesn't skip whitespace */
  82.       if(from.eof()) errfunc(fname); /* find end-of-line */
  83.   }
  84. }
  85.  
  86.  
  87. /* 
  88. -*++ matrix::matrix(): read from a "standard" matrix file
  89. ** 
  90. ** (*++ history: 
  91. **      4 Dec 87    Bruce Eckel    Creation date
  92. ** ++*)
  93. ** 
  94. ** (*++ detailed: 
  95. ** ++*)
  96. */
  97.  
  98.   matrix::matrix(char * initfile) /* read from "standard" matrix file */
  99. {
  100.   void ch_find(istream& from, char ch, void (*errfunc)(char *), char * fname);
  101.   void nonstandard(char *);
  102.   filebuf f1;
  103.   p = new matrep;
  104.   if (f1.open(initfile,input) == 0)
  105.     error2("cannot open matrix initializer file",initfile);
  106.   istream from(&f1);
  107.  
  108.     /* Parse file initialization header  */
  109.  
  110.   ch_find(from,'r',nonstandard,initfile); /* find the 'r'... */
  111.   ch_find(from,'=',nonstandard, initfile); /* ...followed by an `=` */
  112.   from >> p->m_rows;        /* get row size */
  113.   ch_find(from,'c',nonstandard, initfile); /* find the 'p'... */
  114.   ch_find(from,'=',nonstandard, initfile); /* ...followed by an '=' */
  115.   from >> p->m_cols;        /* get column size */
  116.   ch_find(from,':',nonstandard, initfile); /* data follows ::: */
  117.   ch_find(from,':',nonstandard, initfile); 
  118.   ch_find(from,':',nonstandard, initfile); 
  119.   
  120.   p->m_vec = new _Vector*[rows()];
  121.   for (int i = 0; i < rows(); i++)
  122.     p->m_vec[i] = new _Vector(cols());
  123.   p->n = 1;            /* we just created it: only 1 reference */
  124.   for (int row = 0; row < rows(); row++)
  125.     for(int col = 0; col < cols(); col++){
  126.       from >> (*this)[row][col];
  127.       if(from.bad()) error2("problem with matrix initializer file",initfile);
  128.     }
  129. }
  130.  
  131.  
  132. /* 
  133. -*++ nonstandard(): an error message when a "non-standard" matrix file is found
  134. ** 
  135. ** (*++ history: 
  136. **      4 Dec 87    Bruce Eckel    Creation date
  137. ** ++*)
  138. ** 
  139. ** (*++ detailed: 
  140. ** ++*)
  141. */
  142.  
  143. void nonstandard(char * fname){
  144.   cout << fname << " is a `non-standard' file.\n";
  145.   cout << "A `standard' matrix file must start with\n";
  146.   cout << "the dimensions of the matrix, i.e.:\n";
  147.   cout << "rows=12 columns=14\n";
  148.   cout << "or abbreviated:\n";
  149.   cout << "r=12 c=14\n";
  150.   cout << "comments follow `#' signs to end of line\n";
  151.   cout << "data follows :::\n";
  152.   exit(1);
  153. }
  154.  
  155.  
  156. /* 
  157. -*++ matrix::matrix(): create a special square matrix
  158. ** 
  159. ** (*++ history: 
  160. **      4 Dec 87    Bruce Eckel    Creation date
  161. ** ++*)
  162. ** 
  163. ** (*++ detailed: create a special square matrix.  If the "flag" string starts
  164. ** with: 
  165. **   I: identity matrix
  166. **   R: square matrix with random values
  167. ** ++*)
  168. */
  169.  
  170.   matrix::matrix(char * flag, int dimension) 
  171. {
  172.   double drand48();
  173.   int i,j;
  174.   if (flag[0] != 'I' && flag[0] != 'R') 
  175.     error("to create identity: \"I\"; random: \"R\"");
  176.   p = new matrep;
  177.   p->m_vec = new _Vector*[dimension];
  178.   p->m_rows = dimension;
  179.   p->m_cols = dimension;
  180.   for ( i = 0; i < dimension; i++)
  181.     p->m_vec[i] = new _Vector(dimension);
  182.   p->n = 1;            /* so far, only one reference to this data */
  183.  
  184.   switch(flag[0]) {
  185.   case 'I': 
  186.     for (i=0; i < rows(); i++)
  187.       for ( j = 0; j < cols(); j++)
  188.     (*this)[i][j] = (i == j ? 1 : 0);
  189.     break;
  190.   case 'R':
  191.     for ( i=0; i < rows(); i++)
  192.       for ( j = 0; j < cols(); j++)
  193.     (*this)[i][j] = drand48();
  194.   }
  195. }
  196.  
  197.  
  198. /* 
  199. -*++ matrix::~matrix(): matrix destructor.  Probably needs examining
  200. ** 
  201. ** (*++ history: 
  202. **      4 Dec 87    Bruce Eckel    Creation date
  203. ** ++*)
  204. ** 
  205. ** (*++ detailed: 
  206. ** ++*)
  207. */
  208.  
  209.  matrix::~matrix()   // this is probably wrong.
  210. {
  211.   if (--p->n == 0) {
  212.     delete p->m_vec;        /* delete data */
  213.     delete p;
  214.   }
  215. }
  216.  
  217.  
  218. /* 
  219. -*++ matrix::write_standard(): write to a file in "standard" format
  220. ** 
  221. ** (*++ history: 
  222. **      4 Dec 87    Bruce Eckel    Creation date
  223. ** ++*)
  224. ** 
  225. ** (*++ detailed: write this matrix to filename in the 'standard' format
  226. ** ++*)
  227. */
  228.  
  229. void matrix::write_standard(char * filename, char * msg)
  230. {
  231. #ifndef DOS
  232.   char *getwd(char *pathname);
  233. #else DOS
  234.   char *getcwd(char *pathname);
  235. #endif
  236.   char pathname[100];        /* for working directory information */
  237.   filebuf f1;
  238.   if (f1.open(filename,output) == 0)
  239.     error2("cannot open or create matrix output file",filename);
  240.   ostream to(&f1);
  241.   to << "# " << filename << ": matrix file " << " written in `standard' format\n";
  242.   long clock = time(0);
  243.   to << "# " << asctime(localtime(&clock));
  244. #ifndef DOS
  245.   to << "# working directory: " << getwd(pathname) << "\n";
  246. #else DOS
  247.   to << "# working directory: " << getcwd(pathname) << "\n";
  248. #endif
  249.   to << "# " << msg << "\n";
  250.   to << "\nRows=" << rows() << "  Columns=" << cols() << "\n";
  251.   to << ":::\n";
  252.   for (int row = 0; row < rows(); row++){
  253.     for(int col = 0; col < cols(); col++){
  254.       to << form("%6.6g  ",(*this)[row][col]);
  255.       if(to.bad()) error2("problem with matrix output file",filename);
  256.     }
  257.     to << "\n";
  258.   }
  259. }
  260.  
  261.  
  262. /* 
  263. -*++ matrix::operator=(): matrix "equals" using references
  264. ** 
  265. ** (*++ history: 
  266. **      4 Dec 87    Bruce Eckel    Creation date
  267. ** ++*)
  268. ** 
  269. ** (*++ detailed: 
  270. ** ++*)
  271. */
  272.  
  273. matrix &  matrix::operator=(matrix& rval)
  274. {
  275.   rval.p->n++;            /* tell the rval it has another reference to it */
  276.   if(--p->n == 0) {        /* If nobody else is referencing us...*/
  277.     delete p->m_vec;        /* ...nobody else knows to make us go away...*/
  278.     delete p;   
  279.   }
  280.   p = rval.p;     /* ...and we're pointing somewhere else (reference into rval) */
  281.   return *this;
  282. }
  283.  
  284.  
  285. /* 
  286. -*++ matrix::operator+(): add two matrices
  287. ** 
  288. ** (*++ history: 
  289. **      4 Dec 87    Bruce Eckel    Creation date
  290. ** ++*)
  291. ** 
  292. ** (*++ detailed: 
  293. ** ++*)
  294. */
  295.  
  296. matrix & matrix::operator+(matrix& arg)
  297. {
  298.   if((rows() != arg.rows()) || (cols() != arg.cols()))
  299.     error("matrices must be of the same dimension for addition!");
  300.   matrix & sum= *new matrix(rows(),cols());
  301.   for (int i=0; i< rows(); i++)
  302.     for (int j = 0; j < cols(); j++)
  303.       sum[i][j] = (*this)[i][j] + arg[i][j];
  304.   return sum;
  305. }
  306.  
  307.  
  308. /* 
  309. -*++ matrix::operator+(): add a double value to each element
  310. ** 
  311. ** (*++ history: 
  312. **      4 Dec 87    Bruce Eckel    Creation date
  313. ** ++*)
  314. ** 
  315. ** (*++ detailed: 
  316. ** ++*)
  317. */
  318.  
  319. matrix & matrix::operator+(double arg)
  320. {
  321.   matrix & sum= *new matrix(rows(),cols());
  322.   for (int i=0; i< rows(); i++)
  323.     for (int j = 0; j < cols(); j++)
  324.       sum[i][j] = (*this)[i][j] + arg;
  325.   return sum;
  326. }
  327.  
  328.  
  329. /* 
  330. -*++ matrix::operator*(): multiply each element by a double value
  331. ** 
  332. ** (*++ history: 
  333. **      4 Dec 87    Bruce Eckel    Creation date
  334. ** ++*)
  335. ** 
  336. ** (*++ detailed: 
  337. ** ++*)
  338. */
  339.  
  340. matrix & matrix::operator*(double arg)
  341. {
  342.   matrix & sum= *new matrix(rows(),cols());
  343.   for (int i=0; i< rows(); i++)
  344.     for (int j = 0; j < cols(); j++)
  345.       sum[i][j] = (*this)[i][j] * arg;
  346.   return sum;
  347. }
  348.  
  349.  
  350. /* 
  351. -*++ matrix::operator*(): matrix multiply
  352. ** 
  353. ** (*++ history: 
  354. **      4 Dec 87    Bruce Eckel    Creation date
  355. **      9 Dec 87    Bruce Eckel    Changed the test -- I had arg1.rows !=
  356. ** arg2.cols, which was backwards
  357. ** ++*)
  358. ** 
  359. ** (*++ detailed: 
  360. ** ++*)
  361. */
  362.  
  363. matrix & matrix::operator*(matrix& arg)
  364. {
  365.   if(cols() != arg.rows())
  366.     error("# rows of second mat must equal # cols of first for multiply!");
  367.   matrix & result= *new matrix(rows(),arg.cols());
  368.   for(int row = 0; row < rows(); row++)
  369.     for(int col = 0; col < arg.cols(); col++){
  370.       double sum = 0;
  371.       for(int i = 0; i < cols(); i++)
  372.       sum += (*this)[row][i] * arg[i][col];
  373.       result[row][col] = sum;
  374.     };
  375.   return result;
  376. }
  377.  
  378.  
  379. /* 
  380. -*++ matrix::operator-(): subtract one matrix from another
  381. ** 
  382. ** (*++ history: 
  383. **      4 Dec 87    Bruce Eckel    Creation date
  384. ** ++*)
  385. ** 
  386. ** (*++ detailed: 
  387. ** ++*)
  388. */
  389.  
  390. matrix & matrix::operator-(matrix& arg)
  391. {
  392.     if((rows() != arg.rows()) || (cols() != arg.cols()))
  393.     error("matrices must be of the same dimension for subtraction!\n");
  394.     matrix & sum= *new matrix(rows(),cols());
  395.     for (int i=0; i< rows(); i++)
  396.     for (int j = 0; j < cols(); j++)
  397.         sum[i][j] = (*this)[i][j] - arg[i][j];
  398.     return sum;
  399. }
  400.  
  401.  
  402. /* 
  403. -*++ matrix::operator-(): unary minus each element in a matrix
  404. ** 
  405. ** (*++ history: 
  406. **      4 Dec 87    Bruce Eckel    Creation date
  407. ** ++*)
  408. ** 
  409. ** (*++ detailed: 
  410. ** ++*)
  411. */
  412.  
  413. matrix & matrix::operator-()
  414. {
  415.     for (int i=0; i< rows(); i++)
  416.     for (int j = 0; j < cols(); j++)
  417.         (*this)[i][j] *= -1;
  418.     return *this;
  419. }
  420.  
  421.  
  422. /* 
  423. -*++ matrix::operator==(): test for matrix eqivalence
  424. ** 
  425. ** (*++ history: 
  426. **      4 Dec 87    Bruce Eckel    Creation date
  427. ** ++*)
  428. ** 
  429. ** (*++ detailed: 
  430. ** ++*)
  431. */
  432.  
  433. int matrix::operator==(matrix& arg)
  434. {
  435.   if((rows() != arg.rows()) || (cols() != arg.cols()))
  436.     error("matrices must be of the same dimension for equivalence test!\n");
  437.   for (int i=0; i< rows(); i++)
  438.     for (int j = 0; j < cols(); j++)
  439.       if( (*this)[i][j] != arg[i][j])
  440.     return 0;
  441.   return ~0;
  442. }
  443.  
  444.  
  445. /* 
  446. -*++ matrix::operator!=(): test for matrix eqivalence
  447. ** 
  448. ** (*++ history: 
  449. **      4 Dec 87    Bruce Eckel    Creation date
  450. ** ++*)
  451. ** 
  452. ** (*++ detailed: 
  453. ** ++*)
  454. */
  455.  
  456. int matrix::operator!=(matrix& arg)
  457. {
  458.   if((rows() != arg.rows()) || (cols() != arg.cols()))
  459.     error("matrices must be of the same dimension for inequivalence!\n");
  460.   for (int i=0; i< rows(); i++)
  461.     for (int j = 0; j < cols(); j++)
  462.       if( (*this)[i][j] == arg[i][j])
  463.     return 0;
  464.   return ~0;
  465. }
  466.  
  467.  
  468. /* 
  469. -*++ matrix::operator<<: print a matrix to standard out
  470. ** 
  471. ** (*++ history: 
  472. **      4 Dec 87    Bruce Eckel    Creation date
  473. ** ++*)
  474. ** 
  475. ** (*++ detailed: 
  476. ** ++*)
  477. */
  478.  
  479. ostream& operator<<(ostream &s, matrix& mat)
  480. {
  481.   for (int i=0; i< mat.rows(); i++){
  482.     for (int j = 0; j < mat.cols(); j++)
  483.       s <<  form("%6.6g  ",mat[i][j]);
  484.     s << "\n";
  485.   }
  486.   return s;
  487. }
  488.  
  489.  
  490. /* 
  491. -*++ matrix::operator>>(): 
  492. ** 
  493. ** (*++ history: 
  494. **      4 Dec 87    Bruce Eckel    Creation date
  495. ** ++*)
  496. ** 
  497. ** (*++ detailed: 
  498. ** ++*)
  499. */
  500.  
  501. istream& operator>>(istream& s, matrix& mat)
  502. {
  503.   double val = 0;
  504.   for (int row = 0; row < mat.rows(); row++)
  505.     for(int col = 0; col < mat.cols(); col++){
  506.       s >> val;
  507.       if (s.bad()) 
  508.     mat.error("bad matrix initialization input\n");
  509.       if (s.eof()) 
  510.     mat.error("end of matrix initialization input before matrix full\n");
  511.       mat[row][col] = val;
  512.     }
  513.   return s;
  514. }
  515.   
  516.  
  517. /* 
  518. -*++ matrix::file_get(): fills a matrix from a `non-standard' ASCII file  
  519. ** 
  520. ** (*++ history: 
  521. **      4 Dec 87    Bruce Eckel    Creation date
  522. ** ++*)
  523. ** 
  524. ** (*++ detailed: fills a matrix of known size from a `non-standard' file of
  525. ** numbers  
  526. ** ++*)
  527. */
  528.  
  529. void matrix::file_get(char * infile)
  530. {
  531.   filebuf f1;
  532.   if (f1.open(infile,input) == 0)
  533.     error2("cannot open matrix initializer file",infile);
  534.   istream from(&f1);
  535.   for (int row = 0; row < rows(); row++)
  536.     for(int col = 0; col < cols(); col++){
  537.       from >> (*this)[row][col];
  538.       if(from.bad()) error2("problem with matrix initializer file",infile);
  539.     }
  540. }
  541.  
  542.  
  543. /* 
  544. -*++ matrix::file_put(): writes a matrix to a "non-standard" ASCII file
  545. ** 
  546. ** (*++ history: 
  547. **      4 Dec 87    Bruce Eckel    Creation date
  548. ** ++*)
  549. ** 
  550. ** (*++ detailed: 
  551. ** ++*)
  552. */
  553.  
  554. void matrix::file_put(char * outfile)
  555. {
  556.   filebuf f1;
  557.   if (f1.open(outfile,output) == 0)
  558.     error2("cannot open or create matrix output file",outfile);
  559.   ostream to(&f1);
  560.   for (int row = 0; row < rows(); row++){
  561.     for(int col = 0; col < cols(); col++){
  562.       to << (*this)[row][col];
  563.       to << "  ";
  564.       if(to.bad()) error2("problem with matrix output file",outfile);
  565.     }
  566.     to << "\n";
  567.   }
  568. }
  569.  
  570.  
  571. /* 
  572. -*++ operator<<(): output for XY_vec
  573. ** 
  574. ** (*++ history: 
  575. **     11 Dec 87    Bruce Eckel    Creation date
  576. ** ++*)
  577. ** 
  578. ** (*++ detailed: 
  579. ** ++*)
  580. */
  581.  
  582. ostream& operator<<(ostream &s, XY_vec& v) {
  583.     s << " X:       Y:\n";
  584.     for(int i = 0; i < v.size(); i++)
  585.     s << form("%6.6g   %6.6g\n",v.X.Indep(i),v.Y[i]);
  586.     return s;
  587. }
  588.  
  589.  
  590. /* 
  591. -*++ XY_vec::floatvect(): return a pointer to an array of floating-point #s
  592. ** 
  593. ** (*++ history: 
  594. **     14 Dec 87    Bruce Eckel    Creation date
  595. ** ++*)
  596. ** 
  597. ** (*++ detailed: What is this terminated by?
  598. ** ++*)
  599. */
  600.  
  601. float *  XY_vec::floatvect()
  602. {
  603.     float * result = new float[2 * size()];
  604.     for (int i = 0; i < size() * 2; i += 2) {
  605.     result[i] = (float)X.Indep(i/2);
  606.     result[i+1] = (float)Y[i/2];
  607.     }
  608.     return result;
  609. }
  610.  
  611.  
  612. /* 
  613. -*++ ColVec::doublvect(): returns a simple array of double numbers
  614. ** 
  615. ** (*++ history: 
  616. **     15 Jan 88    Bruce Eckel    Creation date
  617. ** ++*)
  618. ** 
  619. ** (*++ detailed: 
  620. ** ++*)
  621. */
  622.  
  623. double * ColVec::doublvect()
  624. {
  625.     double * result = new double[size()];
  626.     for (int i = 0; i < size(); i++)
  627.     result[i] = (*this)[i];
  628.     return result;
  629. }
  630.