home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / vecop.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  14KB  |  598 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. /* vecop.c 1.3 8/18/87 */
  28.  
  29. #include    <stdio.h>
  30. #include    "matrix.h"
  31.  
  32. static    char    rcsid[] = "$Id: vecop.c,v 1.3 1994/01/13 04:29:06 des Exp $";
  33.  
  34.  
  35. /* _in_prod -- inner product of two vectors from i0 downwards */
  36. double    _in_prod(a,b,i0)
  37. VEC    *a,*b;
  38. u_int    i0;
  39. {
  40.     u_int    limit;
  41.     /* Real    *a_v, *b_v; */
  42.     /* register Real    sum; */
  43.  
  44.     if ( a==(VEC *)NULL || b==(VEC *)NULL )
  45.         error(E_NULL,"_in_prod");
  46.     limit = min(a->dim,b->dim);
  47.     if ( i0 > limit )
  48.         error(E_BOUNDS,"_in_prod");
  49.  
  50.     return __ip__(&(a->ve[i0]),&(b->ve[i0]),(int)(limit-i0));
  51.     /*****************************************
  52.     a_v = &(a->ve[i0]);        b_v = &(b->ve[i0]);
  53.     for ( i=i0; i<limit; i++ )
  54.         sum += a_v[i]*b_v[i];
  55.         sum += (*a_v++)*(*b_v++);
  56.  
  57.     return (double)sum;
  58.     ******************************************/
  59. }
  60.  
  61. /* sv_mlt -- scalar-vector multiply -- may be in-situ */
  62. VEC    *sv_mlt(scalar,vector,out)
  63. double    scalar;
  64. VEC    *vector,*out;
  65. {
  66.     /* u_int    dim, i; */
  67.     /* Real    *out_ve, *vec_ve; */
  68.  
  69.     if ( vector==(VEC *)NULL )
  70.         error(E_NULL,"sv_mlt");
  71.     if ( out==(VEC *)NULL || out->dim != vector->dim )
  72.         out = v_resize(out,vector->dim);
  73.     if ( scalar == 0.0 )
  74.         return v_zero(out);
  75.     if ( scalar == 1.0 )
  76.         return v_copy(vector,out);
  77.  
  78.     __smlt__(vector->ve,(double)scalar,out->ve,(int)(vector->dim));
  79.     /**************************************************
  80.     dim = vector->dim;
  81.     out_ve = out->ve;    vec_ve = vector->ve;
  82.     for ( i=0; i<dim; i++ )
  83.         out->ve[i] = scalar*vector->ve[i];
  84.         (*out_ve++) = scalar*(*vec_ve++);
  85.     **************************************************/
  86.     return (out);
  87. }
  88.  
  89. /* v_add -- vector addition -- may be in-situ */
  90. VEC    *v_add(vec1,vec2,out)
  91. VEC    *vec1,*vec2,*out;
  92. {
  93.     u_int    dim;
  94.     /* Real    *out_ve, *vec1_ve, *vec2_ve; */
  95.  
  96.     if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL )
  97.         error(E_NULL,"v_add");
  98.     if ( vec1->dim != vec2->dim )
  99.         error(E_SIZES,"v_add");
  100.     if ( out==(VEC *)NULL || out->dim != vec1->dim )
  101.         out = v_resize(out,vec1->dim);
  102.     dim = vec1->dim;
  103.     __add__(vec1->ve,vec2->ve,out->ve,(int)dim);
  104.     /************************************************************
  105.     out_ve = out->ve;    vec1_ve = vec1->ve;    vec2_ve = vec2->ve;
  106.     for ( i=0; i<dim; i++ )
  107.         out->ve[i] = vec1->ve[i]+vec2->ve[i];
  108.         (*out_ve++) = (*vec1_ve++) + (*vec2_ve++);
  109.     ************************************************************/
  110.  
  111.     return (out);
  112. }
  113.  
  114. /* v_mltadd -- scalar/vector multiplication and addition
  115.         -- out = v1 + scale.v2        */
  116. VEC    *v_mltadd(v1,v2,scale,out)
  117. VEC    *v1,*v2,*out;
  118. double    scale;
  119. {
  120.     /* register u_int    dim, i; */
  121.     /* Real    *out_ve, *v1_ve, *v2_ve; */
  122.  
  123.     if ( v1==(VEC *)NULL || v2==(VEC *)NULL )
  124.         error(E_NULL,"v_mltadd");
  125.     if ( v1->dim != v2->dim )
  126.         error(E_SIZES,"v_mltadd");
  127.     if ( scale == 0.0 )
  128.         return v_copy(v1,out);
  129.     if ( scale == 1.0 )
  130.         return v_add(v1,v2,out);
  131.  
  132.     if ( v2 != out )
  133.     {
  134.         tracecatch(out = v_copy(v1,out),"v_mltadd");
  135.  
  136.         /* dim = v1->dim; */
  137.         __mltadd__(out->ve,v2->ve,scale,(int)(v1->dim));
  138.     }
  139.     else
  140.     {
  141.         tracecatch(out = sv_mlt(scale,v2,out),"v_mltadd");
  142.         out = v_add(v1,out,out);
  143.     }
  144.     /************************************************************
  145.     out_ve = out->ve;    v1_ve = v1->ve;        v2_ve = v2->ve;
  146.     for ( i=0; i < dim ; i++ )
  147.         out->ve[i] = v1->ve[i] + scale*v2->ve[i];
  148.         (*out_ve++) = (*v1_ve++) + scale*(*v2_ve++);
  149.     ************************************************************/
  150.  
  151.     return (out);
  152. }
  153.  
  154. /* v_sub -- vector subtraction -- may be in-situ */
  155. VEC    *v_sub(vec1,vec2,out)
  156. VEC    *vec1,*vec2,*out;
  157. {
  158.     /* u_int    i, dim; */
  159.     /* Real    *out_ve, *vec1_ve, *vec2_ve; */
  160.  
  161.     if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL )
  162.         error(E_NULL,"v_sub");
  163.     if ( vec1->dim != vec2->dim )
  164.         error(E_SIZES,"v_sub");
  165.     if ( out==(VEC *)NULL || out->dim != vec1->dim )
  166.         out = v_resize(out,vec1->dim);
  167.  
  168.     __sub__(vec1->ve,vec2->ve,out->ve,(int)(vec1->dim));
  169.     /************************************************************
  170.     dim = vec1->dim;
  171.     out_ve = out->ve;    vec1_ve = vec1->ve;    vec2_ve = vec2->ve;
  172.     for ( i=0; i<dim; i++ )
  173.         out->ve[i] = vec1->ve[i]-vec2->ve[i];
  174.         (*out_ve++) = (*vec1_ve++) - (*vec2_ve++);
  175.     ************************************************************/
  176.  
  177.     return (out);
  178. }
  179.  
  180. /* v_map -- maps function f over components of x: out[i] = f(x[i])
  181.     -- _v_map sets out[i] = f(x[i],params) */
  182. VEC    *v_map(f,x,out)
  183. double    (*f)();
  184. VEC    *x, *out;
  185. {
  186.     Real    *x_ve, *out_ve;
  187.     int    i, dim;
  188.  
  189.     if ( ! x || ! f )
  190.         error(E_NULL,"v_map");
  191.     if ( ! out || out->dim != x->dim )
  192.         out = v_resize(out,x->dim);
  193.  
  194.     dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  195.     for ( i = 0; i < dim; i++ )
  196.         *out_ve++ = (*f)(*x_ve++);
  197.  
  198.     return out;
  199. }
  200.  
  201. VEC    *_v_map(f,params,x,out)
  202. double    (*f)();
  203. VEC    *x, *out;
  204. void    *params;
  205. {
  206.     Real    *x_ve, *out_ve;
  207.     int    i, dim;
  208.  
  209.     if ( ! x || ! f )
  210.         error(E_NULL,"_v_map");
  211.     if ( ! out || out->dim != x->dim )
  212.         out = v_resize(out,x->dim);
  213.  
  214.     dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  215.     for ( i = 0; i < dim; i++ )
  216.         *out_ve++ = (*f)(params,*x_ve++);
  217.  
  218.     return out;
  219. }
  220.  
  221. /* v_lincomb -- returns sum_i a[i].v[i], a[i] real, v[i] vectors */
  222. VEC    *v_lincomb(n,v,a,out)
  223. int    n;    /* number of a's and v's */
  224. Real    a[];
  225. VEC    *v[], *out;
  226. {
  227.     int    i;
  228.  
  229.     if ( ! a || ! v )
  230.         error(E_NULL,"v_lincomb");
  231.     if ( n <= 0 )
  232.         return VNULL;
  233.  
  234.     for ( i = 1; i < n; i++ )
  235.         if ( out == v[i] )
  236.             error(E_INSITU,"v_lincomb");
  237.  
  238.     out = sv_mlt(a[0],v[0],out);
  239.     for ( i = 1; i < n; i++ )
  240.     {
  241.         if ( ! v[i] )
  242.             error(E_NULL,"v_lincomb");
  243.         if ( v[i]->dim != out->dim )
  244.             error(E_SIZES,"v_lincomb");
  245.         out = v_mltadd(out,v[i],a[i],out);
  246.     }
  247.  
  248.     return out;
  249. }
  250.  
  251.  
  252.  
  253. #ifdef ANSI_C
  254.  
  255. /* v_linlist -- linear combinations taken from a list of arguments;
  256.    calling:
  257.       v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  258.    where vi are vectors (VEC *) and ai are numbers (double)
  259. */
  260. VEC  *v_linlist(VEC *out,VEC *v1,double a1,...)
  261. {
  262.    va_list ap;
  263.    VEC *par;
  264.    double a_par;
  265.  
  266.    if ( ! v1 )
  267.      return VNULL;
  268.    
  269.    va_start(ap, a1);
  270.    out = sv_mlt(a1,v1,out);
  271.    
  272.    while (par = va_arg(ap,VEC *)) {   /* NULL ends the list*/
  273.       a_par = va_arg(ap,double);
  274.       if (a_par == 0.0) continue;
  275.       if ( out == par )        
  276.     error(E_INSITU,"v_linlist");
  277.       if ( out->dim != par->dim )    
  278.     error(E_SIZES,"v_linlist");
  279.  
  280.       if (a_par == 1.0)
  281.     out = v_add(out,par,out);
  282.       else if (a_par == -1.0)
  283.     out = v_sub(out,par,out);
  284.       else
  285.     out = v_mltadd(out,par,a_par,out); 
  286.    } 
  287.    
  288.    va_end(ap);
  289.    return out;
  290. }
  291.  
  292. #elif VARARGS
  293.  
  294.  
  295. /* v_linlist -- linear combinations taken from a list of arguments;
  296.    calling:
  297.       v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  298.    where vi are vectors (VEC *) and ai are numbers (double)
  299. */
  300. VEC  *v_linlist(va_alist) va_dcl
  301. {
  302.    va_list ap;
  303.    VEC *par, *out;
  304.    double a_par;
  305.  
  306.    va_start(ap);
  307.    out = va_arg(ap,VEC *);
  308.    par = va_arg(ap,VEC *);
  309.    if ( ! par ) {
  310.       va_end(ap);
  311.       return VNULL;
  312.    }
  313.    
  314.    a_par = va_arg(ap,double);
  315.    out = sv_mlt(a_par,par,out);
  316.    
  317.    while (par = va_arg(ap,VEC *)) {   /* NULL ends the list*/
  318.       a_par = va_arg(ap,double);
  319.       if (a_par == 0.0) continue;
  320.       if ( out == par )        
  321.     error(E_INSITU,"v_linlist");
  322.       if ( out->dim != par->dim )    
  323.     error(E_SIZES,"v_linlist");
  324.  
  325.       if (a_par == 1.0)
  326.     out = v_add(out,par,out);
  327.       else if (a_par == -1.0)
  328.     out = v_sub(out,par,out);
  329.       else
  330.     out = v_mltadd(out,par,a_par,out); 
  331.    } 
  332.    
  333.    va_end(ap);
  334.    return out;
  335. }
  336.  
  337. #endif
  338.   
  339.  
  340.  
  341.  
  342.  
  343. /* v_star -- computes componentwise (Hadamard) product of x1 and x2
  344.     -- result out is returned */
  345. VEC    *v_star(x1, x2, out)
  346. VEC    *x1, *x2, *out;
  347. {
  348.     int        i;
  349.  
  350.     if ( ! x1 || ! x2 )
  351.     error(E_NULL,"v_star");
  352.     if ( x1->dim != x2->dim )
  353.     error(E_SIZES,"v_star");
  354.     out = v_resize(out,x1->dim);
  355.  
  356.     for ( i = 0; i < x1->dim; i++ )
  357.     out->ve[i] = x1->ve[i] * x2->ve[i];
  358.  
  359.     return out;
  360. }
  361.  
  362. /* v_slash -- computes componentwise ratio of x2 and x1
  363.     -- out[i] = x2[i] / x1[i]
  364.     -- if x1[i] == 0 for some i, then raise E_SING error
  365.     -- result out is returned */
  366. VEC    *v_slash(x1, x2, out)
  367. VEC    *x1, *x2, *out;
  368. {
  369.     int        i;
  370.     Real    tmp;
  371.  
  372.     if ( ! x1 || ! x2 )
  373.     error(E_NULL,"v_slash");
  374.     if ( x1->dim != x2->dim )
  375.     error(E_SIZES,"v_slash");
  376.     out = v_resize(out,x1->dim);
  377.  
  378.     for ( i = 0; i < x1->dim; i++ )
  379.     {
  380.     tmp = x1->ve[i];
  381.     if ( tmp == 0.0 )
  382.         error(E_SING,"v_slash");
  383.     out->ve[i] = x2->ve[i] / tmp;
  384.     }
  385.  
  386.     return out;
  387. }
  388.  
  389. /* v_min -- computes minimum component of x, which is returned
  390.     -- also sets min_idx to the index of this minimum */
  391. double    v_min(x, min_idx)
  392. VEC    *x;
  393. int    *min_idx;
  394. {
  395.     int        i, i_min;
  396.     Real    min_val, tmp;
  397.  
  398.     if ( ! x )
  399.     error(E_NULL,"v_min");
  400.     if ( x->dim <= 0 )
  401.     error(E_SIZES,"v_min");
  402.     i_min = 0;
  403.     min_val = x->ve[0];
  404.     for ( i = 1; i < x->dim; i++ )
  405.     {
  406.     tmp = x->ve[i];
  407.     if ( tmp < min_val )
  408.     {
  409.         min_val = tmp;
  410.         i_min = i;
  411.     }
  412.     }
  413.  
  414.     if ( min_idx != NULL )
  415.     *min_idx = i_min;
  416.     return min_val;
  417. }
  418.  
  419. /* v_max -- computes maximum component of x, which is returned
  420.     -- also sets max_idx to the index of this maximum */
  421. double    v_max(x, max_idx)
  422. VEC    *x;
  423. int    *max_idx;
  424. {
  425.     int        i, i_max;
  426.     Real    max_val, tmp;
  427.  
  428.     if ( ! x )
  429.     error(E_NULL,"v_max");
  430.     if ( x->dim <= 0 )
  431.     error(E_SIZES,"v_max");
  432.     i_max = 0;
  433.     max_val = x->ve[0];
  434.     for ( i = 1; i < x->dim; i++ )
  435.     {
  436.     tmp = x->ve[i];
  437.     if ( tmp > max_val )
  438.     {
  439.         max_val = tmp;
  440.         i_max = i;
  441.     }
  442.     }
  443.  
  444.     if ( max_idx != NULL )
  445.     *max_idx = i_max;
  446.     return max_val;
  447. }
  448.  
  449. #define    MAX_STACK    60
  450.  
  451.  
  452. /* v_sort -- sorts vector x, and generates permutation that gives the order
  453.     of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and
  454.     the permutation is order = [2, 0, 1].
  455.     -- if order is NULL on entry then it is ignored
  456.     -- the sorted vector x is returned */
  457. VEC    *v_sort(x, order)
  458. VEC    *x;
  459. PERM    *order;
  460. {
  461.     Real    *x_ve, tmp, v;
  462.     /* int        *order_pe; */
  463.     int        dim, i, j, l, r, tmp_i;
  464.     int        stack[MAX_STACK], sp;
  465.  
  466.     if ( ! x )
  467.     error(E_NULL,"v_sort");
  468.     if ( order != PNULL && order->size != x->dim )
  469.     order = px_resize(order, x->dim);
  470.  
  471.     x_ve = x->ve;
  472.     dim = x->dim;
  473.     if ( order != PNULL )
  474.     px_ident(order);
  475.  
  476.     if ( dim <= 1 )
  477.     return x;
  478.  
  479.     /* using quicksort algorithm in Sedgewick,
  480.        "Algorithms in C", Ch. 9, pp. 118--122 (1990) */
  481.     sp = 0;
  482.     l = 0;    r = dim-1;    v = x_ve[0];
  483.     for ( ; ; )
  484.     {
  485.     while ( r > l )
  486.     {
  487.         /* "i = partition(x_ve,l,r);" */
  488.         v = x_ve[r];
  489.         i = l-1;
  490.         j = r;
  491.         for ( ; ; )
  492.         {
  493.         while ( x_ve[++i] < v )
  494.             ;
  495.         while ( x_ve[--j] > v )
  496.             ;
  497.         if ( i >= j )    break;
  498.         
  499.         tmp = x_ve[i];
  500.         x_ve[i] = x_ve[j];
  501.         x_ve[j] = tmp;
  502.         if ( order != PNULL )
  503.         {
  504.             tmp_i = order->pe[i];
  505.             order->pe[i] = order->pe[j];
  506.             order->pe[j] = tmp_i;
  507.         }
  508.         }
  509.         tmp = x_ve[i];
  510.         x_ve[i] = x_ve[r];
  511.         x_ve[r] = tmp;
  512.         if ( order != PNULL )
  513.         {
  514.         tmp_i = order->pe[i];
  515.         order->pe[i] = order->pe[r];
  516.         order->pe[r] = tmp_i;
  517.         }
  518.  
  519.         if ( i-l > r-i )
  520.         {   stack[sp++] = l;   stack[sp++] = i-1;   l = i+1;   }
  521.         else
  522.         {   stack[sp++] = i+1;   stack[sp++] = r;   r = i-1;   }
  523.     }
  524.  
  525.     /* recursion elimination */
  526.     if ( sp == 0 )
  527.         break;
  528.     r = stack[--sp];
  529.     l = stack[--sp];
  530.     }
  531.  
  532.     return x;
  533. }
  534.  
  535. /* v_sum -- returns sum of entries of a vector */
  536. double    v_sum(x)
  537. VEC    *x;
  538. {
  539.     int        i;
  540.     Real    sum;
  541.  
  542.     if ( ! x )
  543.     error(E_NULL,"v_sum");
  544.  
  545.     sum = 0.0;
  546.     for ( i = 0; i < x->dim; i++ )
  547.     sum += x->ve[i];
  548.  
  549.     return sum;
  550. }
  551.  
  552. /* v_conv -- computes convolution product of two vectors */
  553. VEC    *v_conv(x1, x2, out)
  554. VEC    *x1, *x2, *out;
  555. {
  556.     int        i;
  557.  
  558.     if ( ! x1 || ! x2 )
  559.     error(E_NULL,"v_conv");
  560.     if ( x1 == out || x2 == out )
  561.     error(E_INSITU,"v_conv");
  562.     if ( x1->dim == 0 || x2->dim == 0 )
  563.     return out = v_resize(out,0);
  564.  
  565.     out = v_resize(out,x1->dim + x2->dim - 1);
  566.     v_zero(out);
  567.     for ( i = 0; i < x1->dim; i++ )
  568.     __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim);
  569.  
  570.     return out;
  571. }
  572.  
  573. /* v_pconv -- computes a periodic convolution product
  574.     -- the period is the dimension of x2 */
  575. VEC    *v_pconv(x1, x2, out)
  576. VEC    *x1, *x2, *out;
  577. {
  578.     int        i;
  579.  
  580.     if ( ! x1 || ! x2 )
  581.     error(E_NULL,"v_pconv");
  582.     if ( x1 == out || x2 == out )
  583.     error(E_INSITU,"v_pconv");
  584.     out = v_resize(out,x2->dim);
  585.     if ( x2->dim == 0 )
  586.     return out;
  587.  
  588.     v_zero(out);
  589.     for ( i = 0; i < x1->dim; i++ )
  590.     {
  591.     __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim - i);
  592.     if ( i > 0 )
  593.         __mltadd__(out->ve,&(x2->ve[x2->dim - i]),x1->ve[i],i);
  594.     }
  595.  
  596.     return out;
  597. }
  598.