home *** CD-ROM | disk | FTP | other *** search
/ Power Programming / powerprogramming1994.iso / progtool / microcrn / issue_40.arc / DAIMS.ARC / MATRIX2.CPP < prev    next >
Text File  |  1988-02-10  |  13KB  |  501 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.  
  9.  
  10. /* 
  11. -*++ matrix::set_row(): sets a row to a value
  12. ** 
  13. ** (*++ history: 
  14. **      4 Dec 87    Bruce Eckel    Creation date
  15. ** ++*)
  16. ** 
  17. ** (*++ detailed: 
  18. ** ++*)
  19. */
  20.  
  21. void matrix::set_row(int row, double value)
  22. {
  23.   for(int i = 0; i < cols(); i++)
  24.     (*this)[row][i] = value;
  25. }
  26.  
  27.  
  28. /* 
  29. -*++ matrix::set_column(): sets a column to a value
  30. ** 
  31. ** (*++ history: 
  32. **      4 Dec 87    Bruce Eckel    Creation date
  33. ** ++*)
  34. ** 
  35. ** (*++ detailed: 
  36. ** ++*)
  37. */
  38.  
  39. void matrix::set_column(int column, double value)
  40. {
  41.   for(int i = 0; i < cols(); i++)
  42.     (*this)[i][column] = value;
  43. }
  44.  
  45.  
  46. /* 
  47. -*++ matrix::switch_rows(): switches two rows of a matrix
  48. ** 
  49. ** (*++ history: 
  50. **      4 Dec 87    Bruce Eckel    Creation date
  51. ** Could be optimized.
  52. ** ++*)
  53. ** 
  54. ** (*++ detailed: 
  55. ** ++*)
  56. */
  57.  
  58. void matrix::switch_rows(int row1, int row2)
  59. {
  60.   matrix temp(1,cols());
  61.   for(int col = 0; col < cols(); col++)
  62.     temp[0][col] = (*this)[row1][col]; /* temporarily store row 1 */
  63.   for(col = 0; col < cols(); col++)
  64.     (*this)[row1][col] = (*this)[row2][col]; /* move row2 to row1 */
  65.   for(col = 0; col < cols(); col++)
  66.     (*this)[row2][col] = temp[0][col]; /* move temp to row2 */
  67. }
  68.  
  69.  
  70. /* 
  71. -*++ matrix::switch_columns(): switches two columns of a matrix
  72. ** 
  73. ** (*++ history: 
  74. **      4 Dec 87    Bruce Eckel    Creation date
  75. ** ++*)
  76. ** 
  77. ** (*++ detailed: 
  78. ** ++*)
  79. */
  80.  
  81. void matrix::switch_columns(int col1, int col2)
  82. {
  83.   matrix temp(rows(),1);
  84.   for(int row = 0; row < rows(); row++)
  85.     temp[row][0] = (*this)[row][col1]; /* temporarily store col 1 */
  86.   for(row = 0; row < rows(); row++)
  87.     (*this)[row][col1] = (*this)[row][col2]; /* move col2 to col1 */
  88.   for(row = 0; row < rows(); row++)
  89.     (*this)[row][col2] = temp[row][0]; /* move temp to col2 */
  90. }
  91.  
  92.  
  93. /* 
  94. -*++ matrix::row_multiply_add(): adds a multiple of one row to another 
  95. ** 
  96. ** (*++ history: 
  97. **      4 Dec 87    Bruce Eckel    Creation date
  98. ** ++*)
  99. ** 
  100. ** (*++ detailed: 
  101. ** ++*)
  102. */
  103.  
  104. void matrix::row_multiply_add(double multiplier, int reference_row, int changing_row)
  105. {
  106.   for(int col = 0; col < cols(); col++)
  107.     (*this)[changing_row][col] += (*this)[reference_row][col] * multiplier;
  108. }
  109.  
  110.  
  111. /* 
  112. -*++ matrix::col_multiply_add(): adds a multiple of one col to another 
  113. ** 
  114. ** (*++ history: 
  115. **      4 Dec 87    Bruce Eckel    Creation date
  116. ** ++*)
  117. ** 
  118. ** (*++ detailed: 
  119. ** ++*)
  120. */
  121.  
  122. void matrix::col_multiply_add(double multiplier, int reference_col, int changing_col)
  123. {
  124.   for(int row = 0; row < rows(); row++)
  125.     (*this)[row][changing_col] += (*this)[row][reference_col] * multiplier;
  126. }
  127.  
  128.  
  129. /* 
  130. -*++ matrix::min(): returns the minimum element in the matrix
  131. ** 
  132. ** (*++ history: 
  133. **      4 Dec 87    Bruce Eckel    Creation date
  134. ** ++*)
  135. ** 
  136. ** (*++ detailed: 
  137. ** ++*)
  138. */
  139.  
  140. double matrix::min()
  141. {
  142.   double temp;
  143.   if(rows() <= 0 || cols() <= 0) error("bad matrix size for min()");
  144.   double minimum = (*this)[0][0];
  145.   for (int row = 0; row < rows(); row++)
  146.     for(int col = 0; col < cols(); col++)
  147.       if ((temp = (*this)[row][col]) < minimum)
  148.     minimum = temp;
  149.   return minimum;
  150. }
  151.  
  152.  
  153. /* 
  154. -*++ matrix::max(): returns the maximum element in the matrix
  155. ** 
  156. ** (*++ history: 
  157. **      4 Dec 87    Bruce Eckel    Creation date
  158. ** ++*)
  159. ** 
  160. ** (*++ detailed: 
  161. ** ++*)
  162. */
  163.  
  164. double matrix::max()
  165. {
  166.   double temp;
  167.   if(rows() <= 0 || cols() <= 0) error("bad matrix size for max()");
  168.   double maximum = (*this)[0][0];
  169.   for (int row = 0; row < rows(); row++)
  170.     for(int col = 0; col < cols(); col++){
  171.       if ((temp = (*this)[row][col]) > maximum)
  172.     maximum = temp;
  173.     }
  174.   return maximum;
  175. }
  176.  
  177.  
  178. /* 
  179. -*++ matrix::mean(): returns the mean of all the matrix elements
  180. ** 
  181. ** (*++ history: 
  182. **      4 Dec 87    Bruce Eckel    Creation date
  183. ** ++*)
  184. ** 
  185. ** (*++ detailed: 
  186. ** ++*)
  187. */
  188.  
  189. double matrix::mean()        
  190. {
  191.   double sum = 0;
  192.   for (int row = 0; row < rows(); row++)
  193.     for(int col = 0; col < cols(); col++)
  194.       sum += (*this)[row][col];
  195.   return sum/(row * col);
  196. }
  197.  
  198.  
  199. /* 
  200. -*++ matrix::variance(): returns the statistical variance of all elements
  201. ** 
  202. ** (*++ history: 
  203. **      4 Dec 87    Bruce Eckel    Creation date
  204. ** ++*)
  205. ** 
  206. ** (*++ detailed: 
  207. ** ++*)
  208. */
  209.  
  210. double matrix::variance()    
  211. {
  212.   double s_squared = 0;
  213.   double mn = mean();
  214.   for (int row = 0; row < rows(); row++)
  215.     for(int col = 0; col < cols(); col++){
  216.       double temp = (*this)[row][col] - mn;
  217.       temp *= temp;
  218.       s_squared += temp;
  219.     }
  220.   s_squared /= row * col -1;    /* number of elements minus one */
  221.   return s_squared;
  222. }
  223.  
  224.  
  225. /* 
  226. -*++ matrix::transpose(): returns the transpose matrix (must be square)
  227. ** 
  228. ** (*++ history: 
  229. **      4 Dec 87    Bruce Eckel    Creation date
  230. ** ++*)
  231. ** 
  232. ** (*++ detailed: 
  233. ** ++*)
  234. */
  235.  
  236. matrix &  matrix::transpose()    
  237. {
  238.   if(rows() != cols()) error("matrix must be square to transpose!\n");
  239.   matrix & trans= *new matrix(rows(),cols());
  240.   for (int row = 0; row < rows(); row++)
  241.     for(int col = 0; col < cols(); col++)
  242.       trans[col][row] = (*this)[row][col];
  243.   return trans;
  244. }
  245.  
  246.  
  247. /* 
  248. -*++ matrix::scale(): scale a matrix (used in LU decomposition)
  249. ** 
  250. ** (*++ history: 
  251. **      4 Dec 87    Bruce Eckel    Creation date
  252. ** ++*)
  253. ** 
  254. ** (*++ detailed: 
  255. ** ++*)
  256. */
  257.  
  258. matrix & matrix::scale()    
  259. {
  260.   double temp;
  261.   if(rows() <= 0 || cols() <= 0) error("bad matrix size for scale()");
  262.   if(rows() != cols()) error("matrix must be square for scale()");
  263.   matrix & scale_vector= *new matrix(rows());
  264.   for (int col = 0; col < cols(); col++){
  265.     double maximum = 0;
  266.     for(int row = 0; row < rows(); row++)
  267.       if ((temp = fabs((*this)[row][col])) > maximum)
  268.     maximum = temp;        /* find max column magnitude in this row */
  269.     if(maximum == 0) error("singular matrix in scale()");
  270.     scale_vector[col][0] = 1/maximum; /* save the scaling */
  271.   }
  272.   return scale_vector;
  273. }
  274.  
  275.  
  276. /* 
  277. -*++ matrix::bitcopy(): make an image of a matrix
  278. ** 
  279. ** (*++ history: 
  280. **      4 Dec 87    Bruce Eckel    Creation date
  281. ** ++*)
  282. ** 
  283. ** (*++ detailed: 
  284. ** ++*)
  285. */
  286.  
  287. void matrix::bitcopy(matrix& from, matrix& to) 
  288. {
  289.   if(from.rows() != to.rows() || from.cols() != to.cols())
  290.     error("matrices must be the same dimensions for bitcopy()");
  291.   for(int row = 0; row < from.rows(); row++)
  292.     for(int col = 0; col < from.cols(); col++)
  293.       to[row][col] = from[row][col];
  294. }
  295.  
  296.  
  297. /* 
  298. -*++ matrix::lu_decompose(): returns LU decomposition matrix
  299. ** 
  300. ** (*++ history: 
  301. **      4 Dec 87    Bruce Eckel    Creation date
  302. ** ++*)
  303. ** 
  304. ** (*++ detailed: indx is an output vector which records the row permutation
  305. ** effected by the partial pivoting, d is output as +-1 depending on whether
  306. ** the number of row interchanges was even or odd, respectively.  This routine 
  307. ** is used in combination with lu_back_subst to solve linear equations or
  308. ** invert a matrix. 
  309. **     Not much in the way of comments here because I don't really know what
  310. ** I'm doing.  I'm following "Numerical Recipes: The Art of Scientific
  311. ** Computing," by Press and Flannery, Cambridge Press 1986 pp. 31-39.     
  312. ** ++*)
  313. */
  314.  
  315. matrix & matrix::lu_decompose(matrix& indx, int& d ) 
  316. {
  317.   if(rows() != cols()) error("Matrix must be square to L-U decompose!\n");
  318.   d = 1;            /* parity check */
  319.   int row,col,k,col_max; /* counters */
  320.   double dum;            /* from the book -- I don't know significance */
  321.   double sum;
  322.   double maximum;
  323.   matrix & lu_decomp = *new matrix(rows(), cols());
  324.   bitcopy(*this,lu_decomp);    /* direct copy of the original matrix */
  325. //  lu_decomp = *this + 0;    /* makes a copy of this (gets around ref)  */
  326.   matrix scale_vector(lu_decomp.scale()); /* scale the matrix */
  327.   
  328.   for(row = 0; row < rows(); row++){ /* The loop over columns of Crout's method */
  329.     if (row > 0) {
  330.       for (col = 0; col < row; col++) { /* eqn 2.3.12 except for row=col */
  331.     sum = lu_decomp[row][col];
  332.     if(col > 0) {
  333.       for(k = 0; k < col; k++)
  334.         sum -= lu_decomp[row][k] * lu_decomp[k][col];
  335.       lu_decomp[row][col] = sum;
  336.     }
  337.       }
  338.     }
  339.     maximum = 0; /* Initialize for the search for the largest pivot element */
  340.     for(col=row; col < cols(); col++){ /* i=j of eq 2.3.12 & i=j+1..N of 2.3.13 */
  341.       sum = lu_decomp[row][col];
  342.       if(row > 0){
  343.     for(k=0; k < row; k++)
  344.       sum -=  lu_decomp[k][col] * lu_decomp[row][k];
  345.     lu_decomp[row][col] = sum;
  346.       }
  347.       dum = scale_vector[col][0] * fabs(sum);/*figure of merit for pivot*/
  348.       if (dum >= maximum){ /* is it better than the best so far? */
  349.     col_max = col;
  350.     maximum = dum;
  351.       }
  352.     }
  353.     col--; /* fix -- maybe difference between fortran and C for-loops. */
  354.        /* I think maybe Fortran leaves the counter at the high */
  355.        /* value while C increments it one past the high value */
  356.     if(row != col_max) {       /* Do we need to interchange rows? */
  357.       lu_decomp.switch_columns(col_max,row); /* Yes, do so... */
  358.       d *= -1;               /* ... and change the parity of d */
  359.       dum = scale_vector[col_max][0]; /* also interchange */ 
  360.       scale_vector[col_max][0] = scale_vector[col][0]; /* the scale factor */
  361.       scale_vector[row][0] = dum; /* Keffer did it this way, not P&F */
  362.     }
  363.     indx[row][0] = col_max;
  364.     if(row != rows() -1){  /* Now, finally, divide by the pivot element */
  365.       if(lu_decomp[row][row] == 0) lu_decomp[row][row] = TINY;
  366.       /* If the pivot element is zero the matrix is singular (at least to the 
  367.      precision of the algorithm).  For some applications on singular matrices,
  368.      it is desirable to substitute TINY for zero/ */
  369.       dum = 1/lu_decomp[row][row];
  370.       for(col=row+1; col < cols(); col++)
  371.     lu_decomp[row][col] *= dum;
  372.     }
  373.   }
  374.   if(lu_decomp[rows()-1][ cols()-1] == 0) 
  375.     lu_decomp[rows()-1][cols()-1] = TINY;
  376.   return lu_decomp;
  377. }
  378.  
  379.  
  380.  
  381. /* 
  382. -*++ matrix::lu_back_subst(): Solves the set of N linear equations A*X = B
  383. ** 
  384. ** (*++ history: 
  385. **      4 Dec 87    Bruce Eckel    Creation date
  386. ** ++*)
  387. ** 
  388. ** (*++ detailed: Here "this" is the LU-decomposition of the matrix A,
  389. ** determined by the routine lu_decompose(). Indx is input as the permutation
  390. ** vector returned by lu_decompose().  B is input as the right-hand side
  391. ** vector B, and returns with the solution vector X.  This routine takes into
  392. ** account the possibility that B will begin with many zero elements, so it is
  393. ** efficient for use in matrix inversion.   See pp 36-37, Press and Flannery.
  394. ** ++*)
  395. */
  396.  
  397.     void matrix::lu_back_subst(matrix& indx, matrix& b)
  398. {
  399.   if(rows() != cols()) 
  400.     error ("non-square lu decomposition matrix in lu_back_subst()");
  401.   if(rows() != b.rows())
  402.     error("wrong size B vector passed to lu_back_subst()");
  403.   if(rows() != indx.rows())
  404.     error("wrong size indx vector passed to lu_back_subst()");
  405.   int row,col,ll;
  406.   int ii = 0;
  407.   double sum;
  408.   for(col=0;col < cols(); col++){
  409.     ll= (int)indx[col][0];
  410.     sum = b[ll][0];
  411.     b[ll][0] = b[col][0];
  412.     if (ii >= 0)
  413.       for(row = ii; row < col; row++)
  414.     sum -= (*this)[row][col] * b[row][0];
  415.     else if(sum != 0)
  416.       ii = col;
  417.     b[col][0] = sum;
  418.   }
  419.   for(col = cols() -1; col >= 0; col--){
  420.     sum = b[col][0];
  421.     if (col < cols() -1)
  422.       for (row = col + 1; row < rows(); row++)
  423.     sum -= (*this)[row][col] * b[row][0];
  424.     b[col][0] = sum/(*this)[col][col]; /* store a component of the soln vector X*/
  425.   }
  426. }
  427.  
  428.  
  429. /* 
  430. -*++ matrix::copy_column(): copy one column onto another
  431. ** 
  432. ** (*++ history: 
  433. **      4 Dec 87    Bruce Eckel    Creation date
  434. ** ++*)
  435. ** 
  436. ** (*++ detailed: copy the from_col of m to the to_col of this
  437. ** ++*)
  438. */
  439.  
  440. void matrix::copy_column(matrix& m, int from_col, int to_col)
  441. {
  442.   if(rows() != m.rows()) error("number of rows must be equal for copy_column()");
  443.   for(int row=0; row < rows(); row++)
  444.     (*this)[row][to_col] = m[row][from_col];
  445. }
  446.  
  447.  
  448.  
  449. /* 
  450. -*++ matrix::inverse(): returns the inverse matrix of a square matrix
  451. ** 
  452. ** (*++ history: 
  453. **      4 Dec 87    Bruce Eckel    Creation date
  454. ** ++*)
  455. ** 
  456. ** (*++ detailed: 
  457. ** ++*)
  458. */
  459.  
  460. matrix & matrix::inverse()    
  461. {
  462.     if(rows() != cols()) error("matrix must be square for inverse()");
  463.     matrix Y("I",rows());        /* create an identity matrix */
  464.     matrix indx(cols());        /* create the "index vector" */
  465.     matrix B(cols());        /* see pp 38. */
  466.     int d;
  467.     matrix & decomp = *new matrix(lu_decompose(indx,d)); 
  468.                 /* perform the decomposition once */
  469.     for(int col = 0; col < cols(); col++){
  470.     B.copy_column(Y,col,0);
  471.     decomp.lu_back_subst(indx,B);
  472.     Y.copy_column(B,0,col);
  473.     }
  474.     return Y.transpose();
  475. }
  476.  
  477.  
  478. /* 
  479. -*++ matrix::determinant(): returns the determinant of a square matrix
  480. ** 
  481. ** (*++ history: 
  482. **      4 Dec 87    Bruce Eckel    Creation date
  483. ** ++*)
  484. ** 
  485. ** (*++ detailed: 
  486. ** ++*)
  487. */
  488.  
  489. double matrix::determinant()
  490. {
  491.   if(rows() != cols()) error("matrix must be square for determinant()");
  492.   matrix indx(cols());        /* create the "index vector" */
  493.   matrix B(cols());        /* see pp 38. */
  494.   int d;
  495.   matrix decomp(lu_decompose(indx,d)); /* perform the decomposition once */
  496.   double determinant = d;
  497.   for(int i=0; i < cols() ; i++)
  498.     determinant *= decomp[i][i];
  499.   return determinant;
  500. }
  501.