home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / extras.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  11KB  |  501 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. /*
  28.     Memory port routines: MEM_COPY and MEM_ZERO
  29. */
  30.  
  31. /* For BSD 4.[23] environments: using bcopy() and bzero() */
  32.  
  33. #include "machine.h"
  34.  
  35. #ifndef MEM_COPY
  36. void    MEM_COPY(from,to,len)
  37. char    *from, *to;
  38. int    len;
  39. {
  40.     int        i;
  41.  
  42.     if ( from < to )
  43.     {
  44.     for ( i = 0; i < len; i++ )
  45.         *to++ = *from++;
  46.     }
  47.     else
  48.     {
  49.     from += len;    to += len;
  50.     for ( i = 0; i < len; i++ )
  51.         *(--to) = *(--from);
  52.     }
  53. }
  54. #endif
  55.  
  56. #ifndef MEM_ZERO
  57. void    MEM_ZERO(ptr,len)
  58. char    *ptr;
  59. int    len;
  60. {
  61.     int        i;
  62.  
  63.     for ( i = 0; i < len; i++ )
  64.     *(ptr++) = '\0';
  65. }
  66. #endif
  67.  
  68. /*
  69.     This file contains versions of something approximating the well-known
  70.     BLAS routines in C, suitable for Meschach (hence the `m').
  71.     These are "vanilla" implementations, at least with some consideration
  72.     of the effects of caching and paging, and maybe some loop unrolling
  73.     for register-rich machines
  74. */
  75.  
  76. /*
  77.     Organisation of matrices: it is assumed that matrices are represented
  78.     by Real **'s. To keep flexibility, there is also an "initial
  79.     column" parameter j0, so that the actual elements used are
  80.         A[0][j0],   A[0][j0+1],   ..., A[0][j0+n-1]
  81.         A[1][j0],   A[1][j0+1],   ..., A[1][j0+n-1]
  82.            ..         ..          ...      ..
  83.         A[m-1][j0], A[m-1][j0+1], ..., A[m-1][j0+n-1]
  84. */
  85.  
  86. static char    rcsid[] = "$Id: extras.c,v 1.3 1994/01/13 05:45:36 des Exp $";
  87.  
  88. #include    <math.h>
  89.  
  90. #define    REGISTER_RICH    1
  91.  
  92. /* mblar-1 routines */
  93.  
  94. /* Mscale -- sets x <- alpha.x */
  95. void    Mscale(len,alpha,x)
  96. int    len;
  97. double    alpha;
  98. Real    *x;
  99. {
  100.     register int    i;
  101.  
  102.     for ( i = 0; i < len; i++ )
  103.     x[i] *= alpha;
  104. }
  105.  
  106. /* Mswap -- swaps x and y */
  107. void    Mswap(len,x,y)
  108. int    len;
  109. Real    *x, *y;
  110. {
  111.     register int    i;
  112.     register Real    tmp;
  113.  
  114.     for ( i = 0; i < len; i++ )
  115.     {
  116.     tmp = x[i];
  117.     x[i] = y[i];
  118.     y[i] = tmp;
  119.     }
  120. }
  121.  
  122. /* Mcopy -- copies x to y */
  123. void    Mcopy(len,x,y)
  124. int    len;
  125. Real    *x, *y;
  126. {
  127.     register int    i;
  128.  
  129.     for ( i = 0; i < len; i++ )
  130.     y[i] = x[i];
  131. }
  132.  
  133. /* Maxpy -- y <- y + alpha.x */
  134. void    Maxpy(len,alpha,x,y)
  135. int    len;
  136. double    alpha;
  137. Real    *x, *y;
  138. {
  139.     register int    i, len4;
  140.  
  141.     /****************************************
  142.     for ( i = 0; i < len; i++ )
  143.     y[i] += alpha*x[i];
  144.     ****************************************/
  145.  
  146. #ifdef REGISTER_RICH
  147.     len4 = len / 4;
  148.     len  = len % 4;
  149.     for ( i = 0; i < len4; i++ )
  150.     {
  151.     y[4*i]   += alpha*x[4*i];
  152.     y[4*i+1] += alpha*x[4*i+1];
  153.     y[4*i+2] += alpha*x[4*i+2];
  154.     y[4*i+3] += alpha*x[4*i+3];
  155.     }
  156.     x += 4*len4;    y += 4*len4;
  157. #endif
  158.     for ( i = 0; i < len; i++ )
  159.     y[i] += alpha*x[i];
  160. }
  161.  
  162. /* Mdot -- returns x'.y */
  163. double    Mdot(len,x,y)
  164. int    len;
  165. Real    *x, *y;
  166. {
  167.     register int    i, len4;
  168.     register Real    sum;
  169.  
  170. #ifndef REGISTER_RICH
  171.     sum = 0.0;
  172. #endif
  173.  
  174. #ifdef REGISTER_RICH
  175.     register Real    sum0, sum1, sum2, sum3;
  176.     
  177.     sum0 = sum1 = sum2 = sum3 = 0.0;
  178.     
  179.     len4 = len / 4;
  180.     len  = len % 4;
  181.     
  182.     for ( i = 0; i < len4; i++ )
  183.     {
  184.     sum0 += x[4*i  ]*y[4*i  ];
  185.     sum1 += x[4*i+1]*y[4*i+1];
  186.     sum2 += x[4*i+2]*y[4*i+2];
  187.     sum3 += x[4*i+3]*y[4*i+3];
  188.     }
  189.     sum = sum0 + sum1 + sum2 + sum3;
  190.     x += 4*len4;    y += 4*len4;
  191. #endif
  192.  
  193.     for ( i = 0; i < len; i++ )
  194.     sum += x[i]*y[i];
  195.  
  196.     return sum;
  197. }
  198.  
  199. #ifndef ABS
  200. #define    ABS(x)    ((x) >= 0 ? (x) : -(x))
  201. #endif
  202.  
  203. /* Mnorminf -- returns ||x||_inf */
  204. double    Mnorminf(len,x)
  205. int    len;
  206. Real    *x;
  207. {
  208.     register int    i;
  209.     register Real    tmp, max_val;
  210.  
  211.     max_val = 0.0;
  212.     for ( i = 0; i < len; i++ )
  213.     {
  214.     tmp = ABS(x[i]);
  215.     if ( max_val < tmp )
  216.         max_val = tmp;
  217.     }
  218.  
  219.     return max_val;
  220. }
  221.  
  222. /* Mnorm1 -- returns ||x||_1 */
  223. double    Mnorm1(len,x)
  224. int    len;
  225. Real    *x;
  226. {
  227.     register int    i;
  228.     register Real    sum;
  229.  
  230.     sum = 0.0;
  231.     for ( i = 0; i < len; i++ )
  232.     sum += ABS(x[i]);
  233.  
  234.     return sum;
  235. }
  236.  
  237. /* Mnorm2 -- returns ||x||_2 */
  238. double    Mnorm2(len,x)
  239. int    len;
  240. Real    *x;
  241. {
  242.     register int    i;
  243.     register Real    norm, invnorm, sum, tmp;
  244.  
  245.     norm = Mnorminf(len,x);
  246.     if ( norm == 0.0 )
  247.     return 0.0;
  248.     invnorm = 1.0/norm;
  249.     sum = 0.0;
  250.     for ( i = 0; i < len; i++ )
  251.     {
  252.     tmp = x[i]*invnorm;
  253.     sum += tmp*tmp;
  254.     }
  255.  
  256.     return sum/invnorm;
  257. }
  258.  
  259. /* mblar-2 routines */
  260.  
  261. /* Mmv -- y <- alpha.A.x + beta.y */
  262. void    Mmv(m,n,alpha,A,j0,x,beta,y)
  263. int    m, n, j0;
  264. double    alpha, beta;
  265. Real    **A, *x, *y;
  266. {
  267.     register int    i, j, m4, n4;
  268.     register Real    sum0, sum1, sum2, sum3, tmp0, tmp1, tmp2, tmp3;
  269.     register Real    *dp0, *dp1, *dp2, *dp3;
  270.  
  271.     /****************************************
  272.     for ( i = 0; i < m; i++ )
  273.     y[i] += alpha*Mdot(n,&(A[i][j0]),x);
  274.     ****************************************/
  275.  
  276.     m4 = n4 = 0;
  277.  
  278. #ifdef REGISTER_RICH
  279.     m4 = m / 4;
  280.     m  = m % 4;
  281.     n4 = n / 4;
  282.     n  = n % 4;
  283.  
  284.     for ( i = 0; i < m4; i++ )
  285.     {
  286.     sum0 = sum1 = sum2 = sum3 = 0.0;
  287.     dp0 = &(A[4*i  ][j0]);
  288.     dp1 = &(A[4*i+1][j0]);
  289.     dp2 = &(A[4*i+2][j0]);
  290.     dp3 = &(A[4*i+3][j0]);
  291.  
  292.     for ( j = 0; j < n4; j++ )
  293.     {
  294.         tmp0 = x[4*j  ];
  295.         tmp1 = x[4*j+1];
  296.         tmp2 = x[4*j+2];
  297.         tmp3 = x[4*j+3];
  298.         sum0 = sum0 + dp0[j]*tmp0 + dp0[j+1]*tmp1 +
  299.         dp0[j+2]*tmp2 + dp0[j+3]*tmp3;
  300.         sum1 = sum1 + dp1[j]*tmp0 + dp1[j+1]*tmp1 +
  301.         dp1[j+2]*tmp2 + dp1[j+3]*tmp3;
  302.         sum2 = sum2 + dp2[j]*tmp0 + dp2[j+1]*tmp1 +
  303.         dp2[j+2]*tmp2 + dp2[j+3]*tmp3;
  304.         sum3 = sum3 + dp3[j]*tmp0 + dp3[j+1]*tmp2 +
  305.         dp3[j+2]*tmp2 + dp3[j+3]*tmp3;
  306.     }
  307.     for ( j = 0; j < n; j++ )
  308.     {
  309.         sum0 += dp0[4*n4+j]*x[4*n4+j];
  310.         sum1 += dp1[4*n4+j]*x[4*n4+j];
  311.         sum2 += dp2[4*n4+j]*x[4*n4+j];
  312.         sum3 += dp3[4*n4+j]*x[4*n4+j];
  313.     }
  314.     y[4*i  ] = beta*y[4*i  ] + alpha*sum0;
  315.     y[4*i+1] = beta*y[4*i+1] + alpha*sum1;
  316.     y[4*i+2] = beta*y[4*i+2] + alpha*sum2;
  317.     y[4*i+3] = beta*y[4*i+3] + alpha*sum3;
  318.     }
  319. #endif
  320.  
  321.     for ( i = 0; i < m; i++ )
  322.     y[4*m4+i] = beta*y[i] + alpha*Mdot(4*n4+n,&(A[4*m4+i][j0]),x);
  323. }
  324.  
  325. /* Mvm -- y <- alpha.A^T.x + beta.y */
  326. void    Mvm(m,n,alpha,A,j0,x,beta,y)
  327. int    m, n, j0;
  328. double    alpha, beta;
  329. Real    **A, *x, *y;
  330. {
  331.     register int    i, j, m4, n2;
  332.     register Real    *Aref;
  333.     register Real     tmp;
  334.  
  335. #ifdef REGISTER_RICH
  336.     register Real    *Aref0, *Aref1;
  337.     register Real    tmp0, tmp1;
  338.     register Real    yval0, yval1, yval2, yval3;
  339. #endif
  340.  
  341.     if ( beta != 1.0 )
  342.     Mscale(m,beta,y);
  343.     /****************************************
  344.     for ( j = 0; j < n; j++ )
  345.     Maxpy(m,alpha*x[j],&(A[j][j0]),y);
  346.     ****************************************/
  347.     m4 = n2 = 0;
  348.  
  349.     m4 = m / 4;
  350.     m  = m % 4;
  351. #ifdef REGISTER_RICH
  352.     n2 = n / 2;
  353.     n  = n % 2;
  354.  
  355.     for ( j = 0; j < n2; j++ )
  356.     {
  357.     tmp0 = alpha*x[2*j];
  358.     tmp1 = alpha*x[2*j+1];
  359.     Aref0 = &(A[2*j  ][j0]);
  360.     Aref1 = &(A[2*j+1][j0]);
  361.     for ( i = 0; i < m4; i++ )
  362.     {
  363.         yval0 = y[4*i  ] + tmp0*Aref0[4*i  ];
  364.         yval1 = y[4*i+1] + tmp0*Aref0[4*i+1];
  365.         yval2 = y[4*i+2] + tmp0*Aref0[4*i+2];
  366.         yval3 = y[4*i+3] + tmp0*Aref0[4*i+3];
  367.         y[4*i  ] = yval0 + tmp1*Aref1[4*i  ];
  368.         y[4*i+1] = yval1 + tmp1*Aref1[4*i+1];
  369.         y[4*i+2] = yval2 + tmp1*Aref1[4*i+2];
  370.         y[4*i+3] = yval3 + tmp1*Aref1[4*i+3];
  371.     }
  372.     y += 4*m4;    Aref0 += 4*m4;    Aref1 += 4*m4;
  373.     for ( i = 0; i < m; i++ )
  374.         y[i] += tmp0*Aref0[i] + tmp1*Aref1[i];
  375.     }
  376. #endif
  377.  
  378.     for ( j = 0; j < n; j++ )
  379.     {
  380.     tmp = alpha*x[2*n2+j];
  381.     Aref = &(A[2*n2+j][j0]);
  382.     for ( i = 0; i < m4; i++ )
  383.     {
  384.         y[4*i  ] += tmp*Aref[4*i  ];
  385.         y[4*i+1] += tmp*Aref[4*i+1];
  386.         y[4*i+2] += tmp*Aref[4*i+2];
  387.         y[4*i+3] += tmp*Aref[4*i+3];
  388.     }
  389.     y += 4*m4;    Aref += 4*m4;
  390.     for ( i = 0; i < m; i++ )
  391.         y[i] += tmp*Aref[i];
  392.     }
  393. }
  394.  
  395. /* Mupdate -- A <- A + alpha.x.y^T */
  396. void    Mupdate(m,n,alpha,x,y,A,j0)
  397. int    m, n, j0;
  398. double    alpha;
  399. Real    **A, *x, *y;
  400. {
  401.     register int    i, j, n4;
  402.     register Real    *Aref;
  403.     register Real     tmp;
  404.  
  405.     /****************************************
  406.     for ( i = 0; i < m; i++ )
  407.     Maxpy(n,alpha*x[i],y,&(A[i][j0]));
  408.     ****************************************/
  409.  
  410.     n4 = n / 4;
  411.     n  = n % 4;
  412.     for ( i = 0; i < m; i++ )
  413.     {
  414.     tmp = alpha*x[i];
  415.     Aref = &(A[i][j0]);
  416.     for ( j = 0; j < n4; j++ )
  417.     {
  418.         Aref[4*j  ] += tmp*y[4*j  ];
  419.         Aref[4*j+1] += tmp*y[4*j+1];
  420.         Aref[4*j+2] += tmp*y[4*j+2];
  421.         Aref[4*j+3] += tmp*y[4*j+3];
  422.     }
  423.     Aref += 4*n4;    y += 4*n4;
  424.     for ( j = 0; j < n; j++ )
  425.         Aref[j] += tmp*y[j];
  426.     }
  427. }
  428.  
  429. /* mblar-3 routines */
  430.  
  431. /* Mmm -- C <- C + alpha.A.B */
  432. void    Mmm(m,n,p,alpha,A,Aj0,B,Bj0,C,Cj0)
  433. int    m, n, p;    /* C is m x n */
  434. double  alpha;
  435. Real    **A, **B, **C;
  436. int    Aj0, Bj0, Cj0;
  437. {
  438.     register int    i, j, k;
  439.     /* register Real    tmp, sum; */
  440.  
  441.     /****************************************
  442.     for ( i = 0; i < m; i++ )
  443.     for ( k = 0; k < p; k++ )
  444.         Maxpy(n,alpha*A[i][Aj0+k],&(B[k][Bj0]),&(C[i][Cj0]));
  445.     ****************************************/
  446.     for ( i = 0; i < m; i++ )
  447.     Mvm(p,n,alpha,&(A[i][Aj0]),B,Bj0,&(C[i][Cj0]));
  448. }
  449.  
  450. /* Mmtrm -- C <- C + alpha.A^T.B */
  451. void    Mmtrm(m,n,p,alpha,A,Aj0,B,Bj0,C,Cj0)
  452. int    m, n, p;    /* C is m x n */
  453. double  alpha;
  454. Real    **A, **B, **C;
  455. int    Aj0, Bj0, Cj0;
  456. {
  457.     register int    i, j, k;
  458.  
  459.     /****************************************
  460.     for ( i = 0; i < m; i++ )
  461.     for ( k = 0; k < p; k++ )
  462.         Maxpy(n,alpha*A[k][Aj0+i],&(B[k][Bj0]),&(C[i][Cj0]));
  463.     ****************************************/
  464.     for ( k = 0; k < p; k++ )
  465.     Mupdate(m,n,alpha,&(A[k][Aj0]),&(B[k][Bj0]),C,Cj0);
  466. }
  467.  
  468. /* Mmmtr -- C <- C + alpha.A.B^T */
  469. void    Mmmtr(m,n,p,alpha,A,Aj0,B,Bj0,C,Cj0)
  470. int    m, n, p;    /* C is m x n */
  471. double  alpha;
  472. Real    **A, **B, **C;
  473. int    Aj0, Bj0, Cj0;
  474. {
  475.     register int    i, j, k;
  476.  
  477.     /****************************************
  478.     for ( i = 0; i < m; i++ )
  479.     for ( j = 0; j < n; j++ )
  480.         C[i][Cj0+j] += alpha*Mdot(p,&(A[i][Aj0]),&(B[j][Bj0]));
  481.     ****************************************/
  482.     for ( i = 0; i < m; i++ )
  483.     Mmv(n,p,alpha,&(A[i][Aj0]),B,Bj0,&(C[i][Cj0]));
  484. }
  485.  
  486. /* Mmtrmtr -- C <- C + alpha.A^T.B^T */
  487. void    Mmtrmtr(m,n,p,alpha,A,Aj0,B,Bj0,C,Cj0)
  488. int    m, n, p;    /* C is m x n */
  489. double  alpha;
  490. Real    **A, **B, **C;
  491. int    Aj0, Bj0, Cj0;
  492. {
  493.     register int    i, j, k;
  494.  
  495.     for ( i = 0; i < m; i++ )
  496.     for ( j = 0; j < n; j++ )
  497.         for ( k = 0; k < p; k++ )
  498.         C[i][Cj0+j] += A[i][Aj0+k]*B[k][Bj0+j];
  499. }
  500.  
  501.