home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / pxop < prev    next >
Encoding:
Text File  |  1994-03-23  |  7.3 KB  |  358 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. /* pxop.c 1.5 12/03/87 */
  28.  
  29.  
  30. #include    <stdio.h>
  31. #include    "matrix.h"
  32.  
  33. static    char    rcsid[] = "$Id: pxop.c,v 1.5 1994/03/23 23:58:50 des Exp $";
  34.  
  35. /**********************************************************************
  36. Note: A permutation is often interpreted as a matrix
  37.         (i.e. a permutation matrix).
  38.     A permutation px represents a permutation matrix P where
  39.         P[i][j] == 1 if and only if px->pe[i] == j
  40. **********************************************************************/
  41.  
  42.  
  43. /* px_inv -- invert permutation -- in situ
  44.     -- taken from ACM Collected Algorithms #250 */
  45. PERM    *px_inv(px,out)
  46. PERM    *px, *out;
  47. {
  48.     int    i, j, k, n, *p;
  49.     
  50.     out = px_copy(px, out);
  51.     n = out->size;
  52.     p = (int *)(out->pe);
  53.     for ( n--; n>=0; n-- )
  54.     {
  55.     i = p[n];
  56.     if ( i < 0 )    p[n] = -1 - i;
  57.     else if ( i != n )
  58.     {
  59.         k = n;
  60.         while (TRUE)
  61.         {
  62.         if ( i < 0 || i >= out->size )
  63.             error(E_BOUNDS,"px_inv");
  64.         j = p[i];    p[i] = -1 - k;
  65.         if ( j == n )
  66.         {    p[n] = i;    break;        }
  67.         k = i;        i = j;
  68.         }
  69.     }
  70.     }
  71.     return out;
  72. }
  73.  
  74. /* px_mlt -- permutation multiplication (composition) */
  75. PERM    *px_mlt(px1,px2,out)
  76. PERM    *px1,*px2,*out;
  77. {
  78.     u_int    i,size;
  79.     
  80.     if ( px1==(PERM *)NULL || px2==(PERM *)NULL )
  81.     error(E_NULL,"px_mlt");
  82.     if ( px1->size != px2->size )
  83.     error(E_SIZES,"px_mlt");
  84.     if ( px1 == out || px2 == out )
  85.     error(E_INSITU,"px_mlt");
  86.     if ( out==(PERM *)NULL || out->size < px1->size )
  87.     out = px_resize(out,px1->size);
  88.     
  89.     size = px1->size;
  90.     for ( i=0; i<size; i++ )
  91.     if ( px2->pe[i] >= size )
  92.         error(E_BOUNDS,"px_mlt");
  93.     else
  94.         out->pe[i] = px1->pe[px2->pe[i]];
  95.     
  96.     return out;
  97. }
  98.  
  99. /* px_vec -- permute vector */
  100. VEC    *px_vec(px,vector,out)
  101. PERM    *px;
  102. VEC    *vector,*out;
  103. {
  104.     u_int    old_i, i, size, start;
  105.     Real    tmp;
  106.     
  107.     if ( px==(PERM *)NULL || vector==(VEC *)NULL )
  108.     error(E_NULL,"px_vec");
  109.     if ( px->size > vector->dim )
  110.     error(E_SIZES,"px_vec");
  111.     if ( out==(VEC *)NULL || out->dim < vector->dim )
  112.     out = v_resize(out,vector->dim);
  113.     
  114.     size = px->size;
  115.     if ( size == 0 )
  116.     return v_copy(vector,out);
  117.     if ( out != vector )
  118.     {
  119.     for ( i=0; i<size; i++ )
  120.         if ( px->pe[i] >= size )
  121.         error(E_BOUNDS,"px_vec");
  122.         else
  123.         out->ve[i] = vector->ve[px->pe[i]];
  124.     }
  125.     else
  126.     {    /* in situ algorithm */
  127.     start = 0;
  128.     while ( start < size )
  129.     {
  130.         old_i = start;
  131.         i = px->pe[old_i];
  132.         if ( i >= size )
  133.         {
  134.         start++;
  135.         continue;
  136.         }
  137.         tmp = vector->ve[start];
  138.         while ( TRUE )
  139.         {
  140.         vector->ve[old_i] = vector->ve[i];
  141.         px->pe[old_i] = i+size;
  142.         old_i = i;
  143.         i = px->pe[old_i];
  144.         if ( i >= size )
  145.             break;
  146.         if ( i == start )
  147.         {
  148.             vector->ve[old_i] = tmp;
  149.             px->pe[old_i] = i+size;
  150.             break;
  151.         }
  152.         }
  153.         start++;
  154.     }
  155.  
  156.     for ( i = 0; i < size; i++ )
  157.         if ( px->pe[i] < size )
  158.         error(E_BOUNDS,"px_vec");
  159.         else
  160.         px->pe[i] = px->pe[i]-size;
  161.     }
  162.     
  163.     return out;
  164. }
  165.  
  166. /* pxinv_vec -- apply the inverse of px to x, returning the result in out */
  167. VEC    *pxinv_vec(px,x,out)
  168. PERM    *px;
  169. VEC    *x, *out;
  170. {
  171.     u_int    i, size;
  172.     
  173.     if ( ! px || ! x )
  174.     error(E_NULL,"pxinv_vec");
  175.     if ( px->size > x->dim )
  176.     error(E_SIZES,"pxinv_vec");
  177.     /* if ( x == out )
  178.     error(E_INSITU,"pxinv_vec"); */
  179.     if ( ! out || out->dim < x->dim )
  180.     out = v_resize(out,x->dim);
  181.     
  182.     size = px->size;
  183.     if ( size == 0 )
  184.     return v_copy(x,out);
  185.     if ( out != x )
  186.     {
  187.     for ( i=0; i<size; i++ )
  188.         if ( px->pe[i] >= size )
  189.         error(E_BOUNDS,"pxinv_vec");
  190.         else
  191.         out->ve[px->pe[i]] = x->ve[i];
  192.     }
  193.     else
  194.     {    /* in situ algorithm --- cheat's way out */
  195.     px_inv(px,px);
  196.     px_vec(px,x,out);
  197.     px_inv(px,px);
  198.     }
  199.  
  200.     return out;
  201. }
  202.  
  203.  
  204.  
  205. /* px_transp -- transpose elements of permutation
  206.         -- Really multiplying a permutation by a transposition */
  207. PERM    *px_transp(px,i1,i2)
  208. PERM    *px;        /* permutation to transpose */
  209. u_int    i1,i2;        /* elements to transpose */
  210. {
  211.     u_int    temp;
  212.  
  213.     if ( px==(PERM *)NULL )
  214.         error(E_NULL,"px_transp");
  215.  
  216.     if ( i1 < px->size && i2 < px->size )
  217.     {
  218.         temp = px->pe[i1];
  219.         px->pe[i1] = px->pe[i2];
  220.         px->pe[i2] = temp;
  221.     }
  222.  
  223.     return px;
  224. }
  225.  
  226. /* myqsort -- a cheap implementation of Quicksort on integers
  227.         -- returns number of swaps */
  228. static int myqsort(a,num)
  229. int    *a, num;
  230. {
  231.     int    i, j, tmp, v;
  232.     int    numswaps;
  233.  
  234.     numswaps = 0;
  235.     if ( num <= 1 )
  236.         return 0;
  237.  
  238.     i = 0;    j = num;    v = a[0];
  239.     for ( ; ; )
  240.     {
  241.         while ( a[++i] < v )
  242.             ;
  243.         while ( a[--j] > v )
  244.             ;
  245.         if ( i >= j )    break;
  246.  
  247.         tmp = a[i];
  248.         a[i] = a[j];
  249.         a[j] = tmp;
  250.         numswaps++;
  251.     }
  252.  
  253.     tmp = a[0];
  254.     a[0] = a[j];
  255.     a[j] = tmp;
  256.     if ( j != 0 )
  257.         numswaps++;
  258.  
  259.     numswaps += myqsort(&a[0],j);
  260.     numswaps += myqsort(&a[j+1],num-(j+1));
  261.  
  262.     return numswaps;
  263. }
  264.  
  265.  
  266. /* px_sign -- compute the ``sign'' of a permutation = +/-1 where
  267.         px is the product of an even/odd # transpositions */
  268. int    px_sign(px)
  269. PERM    *px;
  270. {
  271.     int    numtransp;
  272.     PERM    *px2;
  273.  
  274.     if ( px==(PERM *)NULL )
  275.         error(E_NULL,"px_sign");
  276.     px2 = px_copy(px,PNULL);
  277.     numtransp = myqsort(px2->pe,px2->size);
  278.     px_free(px2);
  279.  
  280.     return ( numtransp % 2 ) ? -1 : 1;
  281. }
  282.  
  283.  
  284. /* px_cols -- permute columns of matrix A; out = A.px'
  285.     -- May NOT be in situ */
  286. MAT    *px_cols(px,A,out)
  287. PERM    *px;
  288. MAT    *A, *out;
  289. {
  290.     int    i, j, m, n, px_j;
  291.     Real    **A_me, **out_me;
  292. #ifdef ANSI_C
  293.     MAT    *m_get(int, int);
  294. #else
  295.     extern MAT    *m_get();
  296. #endif
  297.  
  298.     if ( ! A || ! px )
  299.         error(E_NULL,"px_cols");
  300.     if ( px->size != A->n )
  301.         error(E_SIZES,"px_cols");
  302.     if ( A == out )
  303.         error(E_INSITU,"px_cols");
  304.     m = A->m;    n = A->n;
  305.     if ( ! out || out->m != m || out->n != n )
  306.         out = m_get(m,n);
  307.     A_me = A->me;    out_me = out->me;
  308.  
  309.     for ( j = 0; j < n; j++ )
  310.     {
  311.         px_j = px->pe[j];
  312.         if ( px_j >= n )
  313.             error(E_BOUNDS,"px_cols");
  314.         for ( i = 0; i < m; i++ )
  315.             out_me[i][px_j] = A_me[i][j];
  316.     }
  317.  
  318.     return out;
  319. }
  320.  
  321. /* px_rows -- permute columns of matrix A; out = px.A
  322.     -- May NOT be in situ */
  323. MAT    *px_rows(px,A,out)
  324. PERM    *px;
  325. MAT    *A, *out;
  326. {
  327.     int    i, j, m, n, px_i;
  328.     Real    **A_me, **out_me;
  329. #ifdef ANSI_C
  330.     MAT    *m_get(int, int);
  331. #else
  332.     extern MAT    *m_get();
  333. #endif
  334.  
  335.     if ( ! A || ! px )
  336.         error(E_NULL,"px_rows");
  337.     if ( px->size != A->m )
  338.         error(E_SIZES,"px_rows");
  339.     if ( A == out )
  340.         error(E_INSITU,"px_rows");
  341.     m = A->m;    n = A->n;
  342.     if ( ! out || out->m != m || out->n != n )
  343.         out = m_get(m,n);
  344.     A_me = A->me;    out_me = out->me;
  345.  
  346.     for ( i = 0; i < m; i++ )
  347.     {
  348.         px_i = px->pe[i];
  349.         if ( px_i >= m )
  350.             error(E_BOUNDS,"px_rows");
  351.         for ( j = 0; j < n; j++ )
  352.             out_me[i][j] = A_me[px_i][j];
  353.     }
  354.  
  355.     return out;
  356. }
  357.  
  358.