home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / matop < prev    next >
Encoding:
Text File  |  1994-01-13  |  11.7 KB  |  496 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /* matop.c 1.3 11/25/87 */
  28.  
  29.  
  30. #include    <stdio.h>
  31. #include    "matrix.h"
  32.  
  33. static    char    rcsid[] = "$Id: matop.c,v 1.3 1994/01/13 05:30:28 des Exp $";
  34.  
  35.  
  36. /* m_add -- matrix addition -- may be in-situ */
  37. MAT    *m_add(mat1,mat2,out)
  38. MAT    *mat1,*mat2,*out;
  39. {
  40.     u_int    m,n,i;
  41.  
  42.     if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  43.         error(E_NULL,"m_add");
  44.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  45.         error(E_SIZES,"m_add");
  46.     if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  47.         out = m_resize(out,mat1->m,mat1->n);
  48.     m = mat1->m;    n = mat1->n;
  49.     for ( i=0; i<m; i++ )
  50.     {
  51.         __add__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  52.         /**************************************************
  53.         for ( j=0; j<n; j++ )
  54.             out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
  55.         **************************************************/
  56.     }
  57.  
  58.     return (out);
  59. }
  60.  
  61. /* m_sub -- matrix subtraction -- may be in-situ */
  62. MAT    *m_sub(mat1,mat2,out)
  63. MAT    *mat1,*mat2,*out;
  64. {
  65.     u_int    m,n,i;
  66.  
  67.     if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  68.         error(E_NULL,"m_sub");
  69.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  70.         error(E_SIZES,"m_sub");
  71.     if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  72.         out = m_resize(out,mat1->m,mat1->n);
  73.     m = mat1->m;    n = mat1->n;
  74.     for ( i=0; i<m; i++ )
  75.     {
  76.         __sub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  77.         /**************************************************
  78.         for ( j=0; j<n; j++ )
  79.             out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
  80.         **************************************************/
  81.     }
  82.  
  83.     return (out);
  84. }
  85.  
  86. /* m_mlt -- matrix-matrix multiplication */
  87. MAT    *m_mlt(A,B,OUT)
  88. MAT    *A,*B,*OUT;
  89. {
  90.     u_int    i, /* j, */ k, m, n, p;
  91.     Real    **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
  92.  
  93.     if ( A==(MAT *)NULL || B==(MAT *)NULL )
  94.         error(E_NULL,"m_mlt");
  95.     if ( A->n != B->m )
  96.         error(E_SIZES,"m_mlt");
  97.     if ( A == OUT || B == OUT )
  98.         error(E_INSITU,"m_mlt");
  99.     m = A->m;    n = A->n;    p = B->n;
  100.     A_v = A->me;        B_v = B->me;
  101.  
  102.     if ( OUT==(MAT *)NULL || OUT->m != A->m || OUT->n != B->n )
  103.         OUT = m_resize(OUT,A->m,B->n);
  104.  
  105. /****************************************************************
  106.     for ( i=0; i<m; i++ )
  107.         for  ( j=0; j<p; j++ )
  108.         {
  109.             sum = 0.0;
  110.             for ( k=0; k<n; k++ )
  111.                 sum += A_v[i][k]*B_v[k][j];
  112.             OUT->me[i][j] = sum;
  113.         }
  114. ****************************************************************/
  115.     m_zero(OUT);
  116.     for ( i=0; i<m; i++ )
  117.         for ( k=0; k<n; k++ )
  118.         {
  119.             if ( A_v[i][k] != 0.0 )
  120.                 __mltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p);
  121.             /**************************************************
  122.             B_row = B_v[k];    OUT_row = OUT->me[i];
  123.             for ( j=0; j<p; j++ )
  124.             (*OUT_row++) += tmp*(*B_row++);
  125.             **************************************************/
  126.         }
  127.  
  128.     return OUT;
  129. }
  130.  
  131. /* mmtr_mlt -- matrix-matrix transposed multiplication
  132.     -- A.B^T is returned, and stored in OUT */
  133. MAT    *mmtr_mlt(A,B,OUT)
  134. MAT    *A, *B, *OUT;
  135. {
  136.     int    i, j, limit;
  137.     /* Real    *A_row, *B_row, sum; */
  138.  
  139.     if ( ! A || ! B )
  140.         error(E_NULL,"mmtr_mlt");
  141.     if ( A == OUT || B == OUT )
  142.         error(E_INSITU,"mmtr_mlt");
  143.     if ( A->n != B->n )
  144.         error(E_SIZES,"mmtr_mlt");
  145.     if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
  146.         OUT = m_resize(OUT,A->m,B->m);
  147.  
  148.     limit = A->n;
  149.     for ( i = 0; i < A->m; i++ )
  150.         for ( j = 0; j < B->m; j++ )
  151.         {
  152.             OUT->me[i][j] = __ip__(A->me[i],B->me[j],(int)limit);
  153.             /**************************************************
  154.             sum = 0.0;
  155.             A_row = A->me[i];
  156.             B_row = B->me[j];
  157.             for ( k = 0; k < limit; k++ )
  158.             sum += (*A_row++)*(*B_row++);
  159.             OUT->me[i][j] = sum;
  160.             **************************************************/
  161.         }
  162.  
  163.     return OUT;
  164. }
  165.  
  166. /* mtrm_mlt -- matrix transposed-matrix multiplication
  167.     -- A^T.B is returned, result stored in OUT */
  168. MAT    *mtrm_mlt(A,B,OUT)
  169. MAT    *A, *B, *OUT;
  170. {
  171.     int    i, k, limit;
  172.     /* Real    *B_row, *OUT_row, multiplier; */
  173.  
  174.     if ( ! A || ! B )
  175.         error(E_NULL,"mmtr_mlt");
  176.     if ( A == OUT || B == OUT )
  177.         error(E_INSITU,"mtrm_mlt");
  178.     if ( A->m != B->m )
  179.         error(E_SIZES,"mmtr_mlt");
  180.     if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
  181.         OUT = m_resize(OUT,A->n,B->n);
  182.  
  183.     limit = B->n;
  184.     m_zero(OUT);
  185.     for ( k = 0; k < A->m; k++ )
  186.         for ( i = 0; i < A->n; i++ )
  187.         {
  188.             if ( A->me[k][i] != 0.0 )
  189.             __mltadd__(OUT->me[i],B->me[k],A->me[k][i],(int)limit);
  190.             /**************************************************
  191.             multiplier = A->me[k][i];
  192.             OUT_row = OUT->me[i];
  193.             B_row   = B->me[k];
  194.             for ( j = 0; j < limit; j++ )
  195.             *(OUT_row++) += multiplier*(*B_row++);
  196.             **************************************************/
  197.         }
  198.  
  199.     return OUT;
  200. }
  201.  
  202. /* mv_mlt -- matrix-vector multiplication 
  203.         -- Note: b is treated as a column vector */
  204. VEC    *mv_mlt(A,b,out)
  205. MAT    *A;
  206. VEC    *b,*out;
  207. {
  208.     u_int    i, m, n;
  209.     Real    **A_v, *b_v /*, *A_row */;
  210.     /* register Real    sum; */
  211.  
  212.     if ( A==(MAT *)NULL || b==(VEC *)NULL )
  213.         error(E_NULL,"mv_mlt");
  214.     if ( A->n != b->dim )
  215.         error(E_SIZES,"mv_mlt");
  216.     if ( b == out )
  217.         error(E_INSITU,"mv_mlt");
  218.     if ( out == (VEC *)NULL || out->dim != A->m )
  219.         out = v_resize(out,A->m);
  220.  
  221.     m = A->m;        n = A->n;
  222.     A_v = A->me;        b_v = b->ve;
  223.     for ( i=0; i<m; i++ )
  224.     {
  225.         /* for ( j=0; j<n; j++ )
  226.             sum += A_v[i][j]*b_v[j]; */
  227.         out->ve[i] = __ip__(A_v[i],b_v,(int)n);
  228.         /**************************************************
  229.         A_row = A_v[i];        b_v = b->ve;
  230.         for ( j=0; j<n; j++ )
  231.             sum += (*A_row++)*(*b_v++);
  232.         out->ve[i] = sum;
  233.         **************************************************/
  234.     }
  235.  
  236.     return out;
  237. }
  238.  
  239. /* sm_mlt -- scalar-matrix multiply -- may be in-situ */
  240. MAT    *sm_mlt(scalar,matrix,out)
  241. double    scalar;
  242. MAT    *matrix,*out;
  243. {
  244.     u_int    m,n,i;
  245.  
  246.     if ( matrix==(MAT *)NULL )
  247.         error(E_NULL,"sm_mlt");
  248.     if ( out==(MAT *)NULL || out->m != matrix->m || out->n != matrix->n )
  249.         out = m_resize(out,matrix->m,matrix->n);
  250.     m = matrix->m;    n = matrix->n;
  251.     for ( i=0; i<m; i++ )
  252.         __smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
  253.         /**************************************************
  254.         for ( j=0; j<n; j++ )
  255.             out->me[i][j] = scalar*matrix->me[i][j];
  256.         **************************************************/
  257.     return (out);
  258. }
  259.  
  260. /* vm_mlt -- vector-matrix multiplication 
  261.         -- Note: b is treated as a row vector */
  262. VEC    *vm_mlt(A,b,out)
  263. MAT    *A;
  264. VEC    *b,*out;
  265. {
  266.     u_int    j,m,n;
  267.     /* Real    sum,**A_v,*b_v; */
  268.  
  269.     if ( A==(MAT *)NULL || b==(VEC *)NULL )
  270.         error(E_NULL,"vm_mlt");
  271.     if ( A->m != b->dim )
  272.         error(E_SIZES,"vm_mlt");
  273.     if ( b == out )
  274.         error(E_INSITU,"vm_mlt");
  275.     if ( out == (VEC *)NULL || out->dim != A->n )
  276.         out = v_resize(out,A->n);
  277.  
  278.     m = A->m;        n = A->n;
  279.  
  280.     v_zero(out);
  281.     for ( j = 0; j < m; j++ )
  282.         if ( b->ve[j] != 0.0 )
  283.             __mltadd__(out->ve,A->me[j],b->ve[j],(int)n);
  284.     /**************************************************
  285.     A_v = A->me;        b_v = b->ve;
  286.     for ( j=0; j<n; j++ )
  287.     {
  288.         sum = 0.0;
  289.         for ( i=0; i<m; i++ )
  290.             sum += b_v[i]*A_v[i][j];
  291.         out->ve[j] = sum;
  292.     }
  293.     **************************************************/
  294.  
  295.     return out;
  296. }
  297.  
  298. /* m_transp -- transpose matrix */
  299. MAT    *m_transp(in,out)
  300. MAT    *in, *out;
  301. {
  302.     int    i, j;
  303.     int    in_situ;
  304.     Real    tmp;
  305.  
  306.     if ( in == (MAT *)NULL )
  307.         error(E_NULL,"m_transp");
  308.     if ( in == out && in->n != in->m )
  309.         error(E_INSITU2,"m_transp");
  310.     in_situ = ( in == out );
  311.     if ( out == (MAT *)NULL || out->m != in->n || out->n != in->m )
  312.         out = m_resize(out,in->n,in->m);
  313.  
  314.     if ( ! in_situ )
  315.         for ( i = 0; i < in->m; i++ )
  316.             for ( j = 0; j < in->n; j++ )
  317.                 out->me[j][i] = in->me[i][j];
  318.     else
  319.         for ( i = 1; i < in->m; i++ )
  320.             for ( j = 0; j < i; j++ )
  321.             {    tmp = in->me[i][j];
  322.                 in->me[i][j] = in->me[j][i];
  323.                 in->me[j][i] = tmp;
  324.             }
  325.  
  326.     return out;
  327. }
  328.  
  329. /* swap_rows -- swaps rows i and j of matrix A upto column lim */
  330. MAT    *swap_rows(A,i,j,lo,hi)
  331. MAT    *A;
  332. int    i, j, lo, hi;
  333. {
  334.     int    k;
  335.     Real    **A_me, tmp;
  336.  
  337.     if ( ! A )
  338.         error(E_NULL,"swap_rows");
  339.     if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
  340.         error(E_SIZES,"swap_rows");
  341.     lo = max(0,lo);
  342.     hi = min(hi,A->n-1);
  343.     A_me = A->me;
  344.  
  345.     for ( k = lo; k <= hi; k++ )
  346.     {
  347.         tmp = A_me[k][i];
  348.         A_me[k][i] = A_me[k][j];
  349.         A_me[k][j] = tmp;
  350.     }
  351.     return A;
  352. }
  353.  
  354. /* swap_cols -- swap columns i and j of matrix A upto row lim */
  355. MAT    *swap_cols(A,i,j,lo,hi)
  356. MAT    *A;
  357. int    i, j, lo, hi;
  358. {
  359.     int    k;
  360.     Real    **A_me, tmp;
  361.  
  362.     if ( ! A )
  363.         error(E_NULL,"swap_cols");
  364.     if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
  365.         error(E_SIZES,"swap_cols");
  366.     lo = max(0,lo);
  367.     hi = min(hi,A->m-1);
  368.     A_me = A->me;
  369.  
  370.     for ( k = lo; k <= hi; k++ )
  371.     {
  372.         tmp = A_me[i][k];
  373.         A_me[i][k] = A_me[j][k];
  374.         A_me[j][k] = tmp;
  375.     }
  376.     return A;
  377. }
  378.  
  379. /* ms_mltadd -- matrix-scalar multiply and add
  380.     -- may be in situ
  381.     -- returns out == A1 + s*A2 */
  382. MAT    *ms_mltadd(A1,A2,s,out)
  383. MAT    *A1, *A2, *out;
  384. double    s;
  385. {
  386.     /* register Real    *A1_e, *A2_e, *out_e; */
  387.     /* register int    j; */
  388.     int    i, m, n;
  389.  
  390.     if ( ! A1 || ! A2 )
  391.         error(E_NULL,"ms_mltadd");
  392.     if ( A1->m != A2->m || A1->n != A2->n )
  393.         error(E_SIZES,"ms_mltadd");
  394.  
  395.     if ( s == 0.0 )
  396.         return m_copy(A1,out);
  397.     if ( s == 1.0 )
  398.         return m_add(A1,A2,out);
  399.  
  400.     tracecatch(out = m_copy(A1,out),"ms_mltadd");
  401.  
  402.     m = A1->m;    n = A1->n;
  403.     for ( i = 0; i < m; i++ )
  404.     {
  405.         __mltadd__(out->me[i],A2->me[i],s,(int)n);
  406.         /**************************************************
  407.         A1_e = A1->me[i];
  408.         A2_e = A2->me[i];
  409.         out_e = out->me[i];
  410.         for ( j = 0; j < n; j++ )
  411.             out_e[j] = A1_e[j] + s*A2_e[j];
  412.         **************************************************/
  413.     }
  414.  
  415.     return out;
  416. }
  417.  
  418. /* mv_mltadd -- matrix-vector multiply and add
  419.     -- may not be in situ
  420.     -- returns out == v1 + alpha*A*v2 */
  421. VEC    *mv_mltadd(v1,v2,A,alpha,out)
  422. VEC    *v1, *v2, *out;
  423. MAT    *A;
  424. double    alpha;
  425. {
  426.     /* register    int    j; */
  427.     int    i, m, n;
  428.     Real    *v2_ve, *out_ve;
  429.  
  430.     if ( ! v1 || ! v2 || ! A )
  431.         error(E_NULL,"mv_mltadd");
  432.     if ( out == v2 )
  433.         error(E_INSITU,"mv_mltadd");
  434.     if ( v1->dim != A->m || v2->dim != A-> n )
  435.         error(E_SIZES,"mv_mltadd");
  436.  
  437.     tracecatch(out = v_copy(v1,out),"mv_mltadd");
  438.  
  439.     v2_ve = v2->ve;    out_ve = out->ve;
  440.     m = A->m;    n = A->n;
  441.  
  442.     if ( alpha == 0.0 )
  443.         return out;
  444.  
  445.     for ( i = 0; i < m; i++ )
  446.     {
  447.         out_ve[i] += alpha*__ip__(A->me[i],v2_ve,(int)n);
  448.         /**************************************************
  449.         A_e = A->me[i];
  450.         sum = 0.0;
  451.         for ( j = 0; j < n; j++ )
  452.             sum += A_e[j]*v2_ve[j];
  453.         out_ve[i] = v1->ve[i] + alpha*sum;
  454.         **************************************************/
  455.     }
  456.  
  457.     return out;
  458. }
  459.  
  460. /* vm_mltadd -- vector-matrix multiply and add
  461.     -- may not be in situ
  462.     -- returns out' == v1' + v2'*A */
  463. VEC    *vm_mltadd(v1,v2,A,alpha,out)
  464. VEC    *v1, *v2, *out;
  465. MAT    *A;
  466. double    alpha;
  467. {
  468.     int    /* i, */ j, m, n;
  469.     Real    tmp, /* *A_e, */ *out_ve;
  470.  
  471.     if ( ! v1 || ! v2 || ! A )
  472.         error(E_NULL,"vm_mltadd");
  473.     if ( v2 == out )
  474.         error(E_INSITU,"vm_mltadd");
  475.     if ( v1->dim != A->n || A->m != v2->dim )
  476.         error(E_SIZES,"vm_mltadd");
  477.  
  478.     tracecatch(out = v_copy(v1,out),"vm_mltadd");
  479.  
  480.     out_ve = out->ve;    m = A->m;    n = A->n;
  481.     for ( j = 0; j < m; j++ )
  482.     {
  483.         tmp = v2->ve[j]*alpha;
  484.         if ( tmp != 0.0 )
  485.             __mltadd__(out_ve,A->me[j],tmp,(int)n);
  486.         /**************************************************
  487.         A_e = A->me[j];
  488.         for ( i = 0; i < n; i++ )
  489.             out_ve[i] += A_e[i]*tmp;
  490.         **************************************************/
  491.     }
  492.  
  493.     return out;
  494. }
  495.  
  496.