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