vecop.c

Meschach 可以解稠密或稀疏线性方程组、计算特征值和特征向量和解最小平方问题

/* v_linlist -- linear combinations taken from a list of arguments;
   calling:
      v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
   where vi are vectors (VEC *) and ai are numbers (double)
*/
VEC  *v_linlist(va_alist) va_dcl
{
   va_list ap;
   VEC *par, *out;
   double a_par;

   va_start(ap);
   out = va_arg(ap,VEC *);
   par = va_arg(ap,VEC *);
   if ( ! par ) {
      va_end(ap);
      return VNULL;
   }
   
   a_par = va_arg(ap,double);
   out = sv_mlt(a_par,par,out);
   
   while (par = va_arg(ap,VEC *)) {   /* NULL ends the list*/
      a_par = va_arg(ap,double);
      if (a_par == 0.0) continue;
      if ( out == par )		
	error(E_INSITU,"v_linlist");
      if ( out->dim != par->dim )	
	error(E_SIZES,"v_linlist");

      if (a_par == 1.0)
	out = v_add(out,par,out);
      else if (a_par == -1.0)
	out = v_sub(out,par,out);
      else
	out = v_mltadd(out,par,a_par,out); 
   } 
   
   va_end(ap);
   return out;
}

#endif
  




/* v_star -- computes componentwise (Hadamard) product of x1 and x2
	-- result out is returned */
#ifndef ANSI_C
VEC	*v_star(x1, x2, out)
VEC	*x1, *x2, *out;
#else
VEC	*v_star(const VEC *x1, const VEC *x2, VEC *out)
#endif
{
    int		i;

    if ( ! x1 || ! x2 )
	error(E_NULL,"v_star");
    if ( x1->dim != x2->dim )
	error(E_SIZES,"v_star");
    out = v_resize(out,x1->dim);

    for ( i = 0; i < x1->dim; i++ )
	out->ve[i] = x1->ve[i] * x2->ve[i];

    return out;
}

/* v_slash -- computes componentwise ratio of x2 and x1
	-- out[i] = x2[i] / x1[i]
	-- if x1[i] == 0 for some i, then raise E_SING error
	-- result out is returned */
#ifndef ANSI_C
VEC	*v_slash(x1, x2, out)
VEC	*x1, *x2, *out;
#else
VEC	*v_slash(const VEC *x1, const VEC *x2, VEC *out)
#endif
{
    int		i;
    Real	tmp;

    if ( ! x1 || ! x2 )
	error(E_NULL,"v_slash");
    if ( x1->dim != x2->dim )
	error(E_SIZES,"v_slash");
    out = v_resize(out,x1->dim);

    for ( i = 0; i < x1->dim; i++ )
    {
	tmp = x1->ve[i];
	if ( tmp == 0.0 )
	    error(E_SING,"v_slash");
	out->ve[i] = x2->ve[i] / tmp;
    }

    return out;
}

/* v_min -- computes minimum component of x, which is returned
	-- also sets min_idx to the index of this minimum */
#ifndef ANSI_C
double	v_min(x, min_idx)
VEC	*x;
int	*min_idx;
#else
double	v_min(const VEC *x, int *min_idx)
#endif
{
    int		i, i_min;
    Real	min_val, tmp;

    if ( ! x )
	error(E_NULL,"v_min");
    if ( x->dim <= 0 )
	error(E_SIZES,"v_min");
    i_min = 0;
    min_val = x->ve[0];
    for ( i = 1; i < x->dim; i++ )
    {
	tmp = x->ve[i];
	if ( tmp < min_val )
	{
	    min_val = tmp;
	    i_min = i;
	}
    }

    if ( min_idx != NULL )
	*min_idx = i_min;
    return min_val;
}

/* v_max -- computes maximum component of x, which is returned
	-- also sets max_idx to the index of this maximum */
#ifndef ANSI_C
double	v_max(x, max_idx)
VEC	*x;
int	*max_idx;
#else
double	v_max(const VEC *x, int *max_idx)
#endif
{
    int		i, i_max;
    Real	max_val, tmp;

    if ( ! x )
	error(E_NULL,"v_max");
    if ( x->dim <= 0 )
	error(E_SIZES,"v_max");
    i_max = 0;
    max_val = x->ve[0];
    for ( i = 1; i < x->dim; i++ )
    {
	tmp = x->ve[i];
	if ( tmp > max_val )
	{
	    max_val = tmp;
	    i_max = i;
	}
    }

    if ( max_idx != NULL )
	*max_idx = i_max;
    return max_val;
}

#define	MAX_STACK	60


/* v_sort -- sorts vector x, and generates permutation that gives the order
	of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and
	the permutation is order = [2, 0, 1].
	-- if order is NULL on entry then it is ignored
	-- the sorted vector x is returned */
#ifndef ANSI_C
VEC	*v_sort(x, order)
VEC	*x;
PERM	*order;
#else
VEC	*v_sort(VEC *x, PERM *order)
#endif
{
    Real	*x_ve, tmp, v;
    /* int		*order_pe; */
    int		dim, i, j, l, r, tmp_i;
    int		stack[MAX_STACK], sp;

    if ( ! x )
	error(E_NULL,"v_sort");
    if ( order != PNULL && order->size != x->dim )
	order = px_resize(order, x->dim);

    x_ve = x->ve;
    dim = x->dim;
    if ( order != PNULL )
	px_ident(order);

    if ( dim <= 1 )
	return x;

    /* using quicksort algorithm in Sedgewick,
       "Algorithms in C", Ch. 9, pp. 118--122 (1990) */
    sp = 0;
    l = 0;	r = dim-1;	v = x_ve[0];
    for ( ; ; )
    {
	while ( r > l )
	{
	    /* "i = partition(x_ve,l,r);" */
	    v = x_ve[r];
	    i = l-1;
	    j = r;
	    for ( ; ; )
	    {
		while ( x_ve[++i] < v )
		    ;
		--j;
		while ( x_ve[j] > v && j != 0 )
		    --j;
		if ( i >= j )	break;
		
		tmp = x_ve[i];
		x_ve[i] = x_ve[j];
		x_ve[j] = tmp;
		if ( order != PNULL )
		{
		    tmp_i = order->pe[i];
		    order->pe[i] = order->pe[j];
		    order->pe[j] = tmp_i;
		}
	    }
	    tmp = x_ve[i];
	    x_ve[i] = x_ve[r];
	    x_ve[r] = tmp;
	    if ( order != PNULL )
	    {
		tmp_i = order->pe[i];
		order->pe[i] = order->pe[r];
		order->pe[r] = tmp_i;
	    }

	    if ( i-l > r-i )
	    {   stack[sp++] = l;   stack[sp++] = i-1;   l = i+1;   }
	    else
	    {   stack[sp++] = i+1;   stack[sp++] = r;   r = i-1;   }
	}

	/* recursion elimination */
	if ( sp == 0 )
	    break;
	r = stack[--sp];
	l = stack[--sp];
    }

    return x;
}

/* v_sum -- returns sum of entries of a vector */
#ifndef ANSI_C
double	v_sum(x)
VEC	*x;
#else
double	v_sum(const VEC *x)
#endif
{
    int		i;
    Real	sum;

    if ( ! x )
	error(E_NULL,"v_sum");

    sum = 0.0;
    for ( i = 0; i < x->dim; i++ )
	sum += x->ve[i];

    return sum;
}

/* v_conv -- computes convolution product of two vectors */
#ifndef ANSI_C
VEC	*v_conv(x1, x2, out)
VEC	*x1, *x2, *out;
#else
VEC	*v_conv(const VEC *x1, const VEC *x2, VEC *out)
#endif
{
    int		i;

    if ( ! x1 || ! x2 )
	error(E_NULL,"v_conv");
    if ( x1 == out || x2 == out )
	error(E_INSITU,"v_conv");
    if ( x1->dim == 0 || x2->dim == 0 )
	return out = v_resize(out,0);

    out = v_resize(out,x1->dim + x2->dim - 1);
    v_zero(out);
    for ( i = 0; i < x1->dim; i++ )
	__mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim);

    return out;
}

/* v_pconv -- computes a periodic convolution product
	-- the period is the dimension of x2 */
#ifndef ANSI_C
VEC	*v_pconv(x1, x2, out)
VEC	*x1, *x2, *out;
#else
VEC	*v_pconv(const VEC *x1, const VEC *x2, VEC *out)
#endif
{
    int		i;

    if ( ! x1 || ! x2 )
	error(E_NULL,"v_pconv");
    if ( x1 == out || x2 == out )
	error(E_INSITU,"v_pconv");
    out = v_resize(out,x2->dim);
    if ( x2->dim == 0 )
	return out;

    v_zero(out);
    for ( i = 0; i < x1->dim; i++ )
    {
	__mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim - i);
	if ( i > 0 )
	    __mltadd__(out->ve,&(x2->ve[x2->dim - i]),x1->ve[i],i);
    }

    return out;
}