home *** CD-ROM | disk | FTP | other *** search
/ Shareware Overload / ShartewareOverload.cdr / progm / elem-c.zip / MATRIX.S < prev    next >
Internet Message Format  |  1985-05-18  |  12KB
  1. From hplabs!hao!seismo!rlgvax!cvl!umcp-cs!fred Wed Dec 31 16:00:00 1969
  2. Relay-Version: version B 2.10.1     9/27/83; site hplabsb.UUCP
  3. Posting-Version: version B 2.10 5/3/83; site umcp-cs.UUCP
  4. Path: hplabsb!hplabs!hao!seismo!rlgvax!cvl!umcp-cs!fred
  5. From: fred@umcp-cs.UUCP
  6. Newsgroups: net.sources
  7. Subject: matrix manipulation subroutine library in C
  8. Message-ID: <1184@umcp-cs.UUCP>
  9. Date: Wed, 27-Jul-83 14:35:25 PST
  10. Article-I.D.: umcp-cs.1184
  11. Posted: Wed Jul 27 14:35:25 1983
  12. Date-Received: Thu, 28-Jul-83 00:29:01 PDT
  13. Organization: Univ. of Maryland, Computer Science Dept.
  14. Lines: 472
  15.  
  16. These might be of general interest. Sorry I don't have a makefile
  17. for it, but there are no complex dependency problems.
  18.  
  19. : Run this shell script with "sh" not "csh"
  20. PATH=:/bin:/usr/bin:/usr/ucb
  21. export PATH
  22. all=FALSE
  23. if [ $1x = -ax ]; then
  24.     all=TRUE
  25. fi
  26. /bin/echo 'Extracting matrix.3l'
  27. sed 's/^X//' <<'//go.sysin dd *' >matrix.3l
  28. X.TH MATRIX 3L 12-Oct-1982
  29. X.SH NAME
  30. m_cofactor,  m_copy, m_determ, m_invert, m_multiply, m_solve,
  31. m_transpose \- matrix manipulation subroutines
  32. X.SH SYNOPSIS
  33. X.B #include <local/mat.h>
  34. X.br
  35. X.B cc "*.c" -lmatrix
  36. X.PP
  37. X.B double m_cofactor(m, i, j);
  38. X.br
  39. X.B struct matrix *m;
  40. X.br
  41. X.B int i, j;
  42. X.PP
  43. X.B struct matrix *m_copy(m);
  44. X.br
  45. X.B struct matrix *m;
  46. X.br
  47. X.PP
  48. X.B double m_determ(m);
  49. X.br
  50. X.B struct matrix *m;
  51. X.PP
  52. X.B struct matrix *m_invert(m);
  53. X.br
  54. X.B struct matrix *m;
  55. X.PP
  56. X.B struct matrix *m_multiply(m1, m2);
  57. X.br
  58. X.B struct matrix *m1;
  59. X.br
  60. X.B struct matrix *m2;
  61. X.PP
  62. X.B struct matrix *m_solve(a, b);
  63. X.br
  64. X.B struct matrix *a;
  65. X.br
  66. X.B struct matrix *b;
  67. X.PP
  68. X.B struct matrix *m_transpose(m);
  69. X.br
  70. X.B struct matrix *m;
  71. X.SH DESCRIPTION
  72. This set of subroutines does several useful operations on matrices which are
  73. stored in a standard format (defined in /usr/include/local/matrix.h).
  74. Space for matrices is allocated and freed by the subroutines
  75. X.I malloc
  76. and
  77. X.I free.
  78. Matrices should be initialized by the macro
  79. X.I m_create,
  80. and can be disposed of by passing a pointer to the matrix to
  81. X.I free.
  82. X.PP
  83. X.I m_cofactor\c
  84. (m, i, j) returns the cofactor, as type
  85. X.I double,
  86. of element (i, j) of the matrix given by its first argument.
  87. That is: the determinant of the matrix obtained by deleting row
  88. X.I i
  89. and column
  90. X.I j
  91. from the matrix
  92. X.I m,
  93. multiplied by
  94. X.br
  95. (-1)**(i+j).
  96. X.PP
  97. X.I m_copy\c
  98. (m) returns a duplicate copy of the matrix
  99. X.I m.
  100. X.PP
  101. X.I m_determ\c
  102. (m) returns the determinant of the matrix m as type
  103. X.I double.
  104. X.PP
  105. X.I m_invert\c
  106. (m) returns a matrix (if it exists) which is the inverse of the matrix
  107. X.I m.
  108. X.PP
  109. X.I m_multiply\c
  110. (m1, m2) returns a matrix which is the result of multiplying matrix
  111. X.I m1
  112. by matrix
  113. X.I m2
  114. by the conventional definition of matrix multiplication. The number of
  115. columns in
  116. X.I m1
  117. must be the same as the number of rows in
  118. X.I m2.
  119. X.PP
  120. X.I m_solve\c
  121. (a, b) returns the matrix which results from computing: (a'\ a)**(-1)\ a'\ b.
  122. This is useful for solving a set of linear equations expressed as matrices:
  123. a\ x\ =\ b.
  124. X.PP
  125. X.I m_transpose\c
  126. (m) returns the matrix which is the transpose of
  127. X.I m.
  128. X.PP
  129. X.I m_create\c
  130. (m, r, c) (which is a macro, not a function) allocates space for,
  131. and initializes the matrix
  132. X.I
  133. m
  134. to have
  135. X.I r
  136. rows and
  137. X.I c
  138. columns.
  139. X.PP
  140. X.I m_v\c
  141. (m, r, c) is a macro used to access element
  142. X.I (r, c)
  143. of the matrix
  144. X.I m.
  145. The elements are always of type double.
  146. X.SH AUTHOR
  147. Fred Blonder <fred@umcp-cs>
  148. X.SH "SEE ALSO"
  149. malloc(3), free(3), printf(3)
  150. X.SH DIAGNOSTICS
  151. Most routines return the constant
  152. X.I M_NULL
  153. to report an error. Some also write messages to the standard output via
  154. X.I printf.
  155. X.SH FILES
  156. X/usr/include/local/matrix.h
  157. X.br
  158. X/usr/lib/libmatrix.a
  159. X.SH BUGS
  160. Error reporting via
  161. X.I printf
  162. should be an option.
  163. //go.sysin dd *
  164. made=TRUE
  165. if [ $made = TRUE ]; then
  166.     /bin/chmod 644 matrix.3l
  167.     /bin/echo -n '    '; /bin/ls -ld matrix.3l
  168. fi
  169. /bin/echo 'Extracting mat.h'
  170. sed 's/^X//' <<'//go.sysin dd *' >mat.h
  171. struct matrix {
  172.     int m_rows, m_cols;
  173.     };
  174.  
  175. struct m_ {
  176.     int m_rows, m_cols;
  177.     double m_value[1];
  178.     };
  179.  
  180. double m_cofactor(), m_determinant();
  181. struct matrix *m_copy(), *m_invert(), *m_transpose(), *m_multiply(), *m_solve();
  182.  
  183. #define    m_v(m, r, c)    ((struct m_ *)m)->m_value[r * (m->m_cols) + c]
  184. #define    M_NULL    ((struct matrix *)0)
  185.  
  186. #define    m_create(m, r, c) {\
  187.             if (((int)(m = (struct matrix *)malloc(sizeof(struct matrix) + (sizeof(double) * r * c)))) == 0) {\
  188.                 printf("Allocation error: %s\n", __FILE__);\
  189.                 exit(1);\
  190.                 }\
  191.             m->m_rows = r;\
  192.             m->m_cols = c;\
  193.             }
  194. //go.sysin dd *
  195. made=TRUE
  196. if [ $made = TRUE ]; then
  197.     /bin/chmod 644 mat.h
  198.     /bin/echo -n '    '; /bin/ls -ld mat.h
  199. fi
  200. /bin/echo 'Extracting m_cofactor.c'
  201. sed 's/^X//' <<'//go.sysin dd *' >m_cofactor.c
  202. static char *sccsid = "@(#)m_cofactor.c    4/5/82 (U of Maryland, FLB)";
  203.  
  204. #include "mat.h"
  205.  
  206. double
  207. m_cofactor(mat, i, j)
  208. register struct matrix *mat;
  209. register int i, j;
  210. {
  211. register struct matrix *result;
  212. register int row, col, o_row = 0, o_col = 0;
  213. double det;
  214.  
  215. m_create(result, mat->m_rows - 1, mat->m_cols - 1);
  216.  
  217. for (row = 0; row < mat->m_rows; row++)
  218.     for (col = 0; col < mat->m_cols; col++)
  219.         if (row != i && col != j)
  220.             m_v(result, o_row++, o_col++) = m_v(mat, row, col);
  221.  
  222. det = m_determinant(result);
  223. free(result);
  224.  
  225. return(((i + j) & 01)? -det: det);
  226. }
  227. //go.sysin dd *
  228. made=TRUE
  229. if [ $made = TRUE ]; then
  230.     /bin/chmod 644 m_cofactor.c
  231.     /bin/echo -n '    '; /bin/ls -ld m_cofactor.c
  232. fi
  233. /bin/echo 'Extracting m_copy.c'
  234. sed 's/^X//' <<'//go.sysin dd *' >m_copy.c
  235. static char *sccsid = "@(#)m_copy.c    4/5/82 (U of Maryland, FLB)";
  236.  
  237. #include "mat.h"
  238.  
  239. struct matrix *
  240. m_copy(mat)
  241. register struct matrix *mat;
  242. {
  243. register struct matrix *result;
  244. register int row, col;
  245.  
  246. m_create(result, mat->m_rows, mat->m_cols);
  247.  
  248. for (row = 0; row < mat->m_rows; row++)
  249.     for (col = 0; col < mat->m_cols; col++)
  250.         m_v(result, row, col) = m_v(mat, row, col);
  251.  
  252. return(result);
  253. }
  254. //go.sysin dd *
  255. made=TRUE
  256. if [ $made = TRUE ]; then
  257.     /bin/chmod 644 m_copy.c
  258.     /bin/echo -n '    '; /bin/ls -ld m_copy.c
  259. fi
  260. /bin/echo 'Extracting m_determ.c'
  261. sed 's/^X//' <<'//go.sysin dd *' >m_determ.c
  262. static char *sccsid = "@(#)m_determinant.c    4/5/82 (U of Maryland, FLB)";
  263.  
  264. #include "mat.h"
  265.  
  266. double
  267. m_determinant(mat)
  268. register struct matrix *mat;
  269. {
  270. register int col;
  271. double det = 0.0;
  272.  
  273. if (mat->m_rows == 1)
  274.     return(m_v(mat, 0, 0));
  275.  
  276. for (col = 0; col < mat->m_cols; col++)
  277.     det +=     m_v(mat, 0, col) * m_cofactor(mat, 0, col);
  278.  
  279. return(det);
  280. }
  281. //go.sysin dd *
  282. made=TRUE
  283. if [ $made = TRUE ]; then
  284.     /bin/chmod 644 m_determ.c
  285.     /bin/echo -n '    '; /bin/ls -ld m_determ.c
  286. fi
  287. /bin/echo 'Extracting m_dump.c'
  288. sed 's/^X//' <<'//go.sysin dd *' >m_dump.c
  289. static char *sccsid = "@(#)m_dump.c    4/6/82 (U of Maryland, FLB)";
  290.  
  291. #include "mat.h"
  292.  
  293. struct matrix *
  294. m_dump(mat)
  295. register struct matrix *mat;
  296. {
  297. register int row, col;
  298.  
  299. printf("%d X %d\n", mat->m_rows, mat->m_cols);
  300.  
  301. for (row = 0; row < mat->m_rows; row++) {
  302.     for (col = 0; col < mat->m_cols; col++)
  303.         printf("%f, ", m_v(mat, row, col));
  304.     putchar('\n');
  305.     }
  306. }
  307. //go.sysin dd *
  308. made=TRUE
  309. if [ $made = TRUE ]; then
  310.     /bin/chmod 644 m_dump.c
  311.     /bin/echo -n '    '; /bin/ls -ld m_dump.c
  312. fi
  313. /bin/echo 'Extracting m_invert.c'
  314. sed 's/^X//' <<'//go.sysin dd *' >m_invert.c
  315. static char *sccsid = "@(#)m_invert.c    4/7/82 (U of Maryland, FLB)";
  316.  
  317. #include "mat.h"
  318.  
  319. struct matrix *
  320. m_invert(mat)
  321. register struct matrix *mat;
  322. {
  323. register struct matrix *result;
  324. register int row, col;
  325. double det;
  326.  
  327. if((det = m_determinant(mat)) == 0.0) {
  328.     printf("m_invert: singular matrix\n");
  329.     return(M_NULL);
  330.     }
  331.  
  332. m_create(result, mat->m_cols, mat->m_rows);
  333.  
  334. for (row = 0; row < mat->m_rows; row++)
  335.     for (col = 0; col < mat->m_cols; col++)
  336.         m_v(result, col, row) = m_cofactor(mat, row, col) / det;
  337.  
  338. return(result);
  339. }
  340. //go.sysin dd *
  341. made=TRUE
  342. if [ $made = TRUE ]; then
  343.     /bin/chmod 644 m_invert.c
  344.     /bin/echo -n '    '; /bin/ls -ld m_invert.c
  345. fi
  346. /bin/echo 'Extracting m_multiply.c'
  347. sed 's/^X//' <<'//go.sysin dd *' >m_multiply.c
  348. static char *sccsid = "@(#)m_multiply.c    4/6/82 (U of Maryland, FLB)";
  349.  
  350. #include "mat.h"
  351.  
  352. struct matrix *
  353. m_multiply(mat1, mat2)
  354. register struct matrix *mat1, *mat2;
  355. {
  356. register struct matrix *result;
  357. register int row, col, ix;
  358. double sum;
  359.  
  360. if (mat1->m_cols != mat2->m_rows) {
  361.     printf("m_multiply: matrices not compatible.\n");
  362.     return(M_NULL);
  363.     }
  364.  
  365. m_create(result, mat1->m_rows, mat2->m_cols);
  366.  
  367. for (row = 0; row < mat1->m_rows; row++)
  368.     for (col = 0; col < mat2->m_cols; col++) {
  369.         sum = 0.0;
  370.         for (ix = 0; ix < mat1->m_cols; ix++)
  371.             sum += m_v(mat1, row, ix) * m_v(mat2, ix, col);
  372.         m_v(result, row, col) = sum;
  373.         }
  374.  
  375. return(result);
  376. }
  377. //go.sysin dd *
  378. made=TRUE
  379. if [ $made = TRUE ]; then
  380.     /bin/chmod 644 m_multiply.c
  381.     /bin/echo -n '    '; /bin/ls -ld m_multiply.c
  382. fi
  383. /bin/echo 'Extracting m_read.c'
  384. sed 's/^X//' <<'//go.sysin dd *' >m_read.c
  385. static char *sccsid = "@(#)m_read.c    4/6/82 (U of Maryland, FLB)";
  386.  
  387. #include <stdio.h>
  388. #include "mat.h"
  389.  
  390. struct matrix *
  391. m_read()
  392. {
  393. register struct matrix *result;
  394. register int row, col;
  395. int rows, cols;
  396.  
  397. scanf("%d%d", &rows, &cols);
  398. m_create(result, rows, cols);
  399.  
  400. for (row = 0; row < rows; row++)
  401.     for (col = 0; col < cols; col++)
  402.         scanf("%lf", &m_v(result, row, col));
  403.  
  404. return(result);
  405. }
  406. //go.sysin dd *
  407. made=TRUE
  408. if [ $made = TRUE ]; then
  409.     /bin/chmod 644 m_read.c
  410.     /bin/echo -n '    '; /bin/ls -ld m_read.c
  411. fi
  412. /bin/echo 'Extracting m_solve.c'
  413. sed 's/^X//' <<'//go.sysin dd *' >m_solve.c
  414. static char *sccsid = "@(#)m_solve.c    4/16/82 (U of Maryland, FLB)";
  415.  
  416. #include "mat.h"
  417.  
  418. struct matrix *
  419. m_solve(a, b)
  420. register struct matrix *a, *b;
  421. {
  422. register struct matrix *a_trans, *t, *t_inv, *t2, *x;
  423.  
  424. if ((a->m_rows) < (a->m_cols)) {
  425.     printf("m_solve: too few equations\n");
  426.     return(M_NULL);
  427.     }
  428.  
  429. if ((a->m_rows) != (b->m_rows)) {
  430.     printf("m_solve: arguments don't match: %d, %d.\n", a->m_rows, b->m_rows);
  431.     return(M_NULL);
  432.     }
  433.  
  434. if (b->m_cols != 1) {
  435.     printf("m_solve: arg2 must have 1 column.\n");
  436.     return(M_NULL);
  437.     }
  438.  
  439. a_trans = m_transpose(a);        /* A' */
  440. t = m_multiply(a_trans, a);        /* A' A */
  441. t_inv = m_invert(t);            /* (A' A)-1 */
  442. free(t);
  443. if (t_inv == M_NULL) {
  444.     printf("m_solve: no solution\n");
  445.     return(M_NULL);
  446.     }
  447. t2 = m_multiply(t_inv, a_trans);    /* (A' A)-1 A' */
  448. free(t_inv);
  449. free(a_trans);
  450. x = m_multiply(t2, b);            /* (A' A)-1 A' B */
  451. free(t2);
  452.  
  453. return(x);
  454. }
  455. //go.sysin dd *
  456. made=TRUE
  457. if [ $made = TRUE ]; then
  458.     /bin/chmod 644 m_solve.c
  459.     /bin/echo -n '    '; /bin/ls -ld m_solve.c
  460. fi
  461. /bin/echo 'Extracting m_transpose.c'
  462. sed 's/^X//' <<'//go.sysin dd *' >m_transpose.c
  463. static char *sccsid = "@(#)m_transpose.c    4/5/82 (U of Maryland, FLB)";
  464.  
  465. #include "mat.h"
  466.  
  467. struct matrix *
  468. m_transpose(mat)
  469. register struct matrix *mat;
  470. {
  471. register struct matrix *result;
  472. register int row, col;
  473.  
  474. m_create(result, mat->m_cols, mat->m_rows);
  475.  
  476. for (row = 0; row < mat->m_rows; row++)
  477.     for (col = 0; col < mat->m_cols; col++)
  478.         m_v(result, col, row) = m_v(mat, row, col);
  479.  
  480. return(result);
  481. }
  482. //go.sysin dd *
  483. made=TRUE
  484. if [ $made = TRUE ]; then
  485.     /bin/chmod 644 m_transpose.c
  486.     /bin/echo -n '    '; /bin/ls -ld m_transpose.c
  487. fi
  488.  
  489.