zmatop.c

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

    if ( in == out && in->n != in->m )
	error(E_INSITU2,"zm_adjoint");
    in_situ = ( in == out );
    if ( out == ZMNULL || out->m != in->n || out->n != in->m )
	out = zm_resize(out,in->n,in->m);
    
    if ( ! in_situ )
    {
	for ( i = 0; i < in->m; i++ )
	    for ( j = 0; j < in->n; j++ )
	    {
		out->me[j][i].re =   in->me[i][j].re;
		out->me[j][i].im = - in->me[i][j].im;
	    }
    }
    else
    {
	for ( i = 0 ; i < in->m; i++ )
	{
	    for ( j = 0; j < i; j++ )
	    {
		tmp.re = in->me[i][j].re;
		tmp.im = in->me[i][j].im;
		in->me[i][j].re =   in->me[j][i].re;
		in->me[i][j].im = - in->me[j][i].im;
		in->me[j][i].re =   tmp.re;
		in->me[j][i].im = - tmp.im;
	    }
	    in->me[i][i].im = - in->me[i][i].im;
	}
    }
    
    return out;
}

/* zswap_rows -- swaps rows i and j of matrix A upto column lim */
ZMAT	*zswap_rows(A,i,j,lo,hi)
ZMAT	*A;
int	i, j, lo, hi;
{
    int	k;
    complex	**A_me, tmp;
    
    if ( ! A )
	error(E_NULL,"swap_rows");
    if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
	error(E_SIZES,"swap_rows");
    lo = max(0,lo);
    hi = min(hi,A->n-1);
    A_me = A->me;
    
    for ( k = lo; k <= hi; k++ )
    {
	tmp = A_me[k][i];
	A_me[k][i] = A_me[k][j];
	A_me[k][j] = tmp;
    }
    return A;
}

/* zswap_cols -- swap columns i and j of matrix A upto row lim */
ZMAT	*zswap_cols(A,i,j,lo,hi)
ZMAT	*A;
int	i, j, lo, hi;
{
    int	k;
    complex	**A_me, tmp;
    
    if ( ! A )
	error(E_NULL,"swap_cols");
    if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
	error(E_SIZES,"swap_cols");
    lo = max(0,lo);
    hi = min(hi,A->m-1);
    A_me = A->me;
    
    for ( k = lo; k <= hi; k++ )
    {
	tmp = A_me[i][k];
	A_me[i][k] = A_me[j][k];
	A_me[j][k] = tmp;
    }
    return A;
}

/* mz_mltadd -- matrix-scalar multiply and add
   -- may be in situ
   -- returns out == A1 + s*A2 */
ZMAT	*mz_mltadd(A1,A2,s,out)
ZMAT	*A1, *A2, *out;
complex	s;
{
    /* register complex	*A1_e, *A2_e, *out_e; */
    /* register int	j; */
    int	i, m, n;
    
    if ( ! A1 || ! A2 )
	error(E_NULL,"mz_mltadd");
    if ( A1->m != A2->m || A1->n != A2->n )
	error(E_SIZES,"mz_mltadd");

    if ( out != A1 && out != A2 )
        out = zm_resize(out,A1->m,A1->n);
    
    if ( s.re == 0.0 && s.im == 0.0 )
	return zm_copy(A1,out);
    if ( s.re == 1.0 && s.im == 0.0 )
	return zm_add(A1,A2,out);
    
    out = zm_copy(A1,out);
    
    m = A1->m;	n = A1->n;
    for ( i = 0; i < m; i++ )
    {
	__zmltadd__(out->me[i],A2->me[i],s,(int)n,Z_NOCONJ);
	/**************************************************
	  A1_e = A1->me[i];
	  A2_e = A2->me[i];
	  out_e = out->me[i];
	  for ( j = 0; j < n; j++ )
	  out_e[j] = A1_e[j] + s*A2_e[j];
	  **************************************************/
    }
    
    return out;
}

/* zmv_mltadd -- matrix-vector multiply and add
   -- may not be in situ
   -- returns out == v1 + alpha*A*v2 */
ZVEC	*zmv_mltadd(v1,v2,A,alpha,out)
ZVEC	*v1, *v2, *out;
ZMAT	*A;
complex	alpha;
{
    /* register	int	j; */
    int	i, m, n;
    complex	tmp, *v2_ve, *out_ve;
    
    if ( ! v1 || ! v2 || ! A )
	error(E_NULL,"zmv_mltadd");
    if ( out == v2 )
	error(E_INSITU,"zmv_mltadd");
    if ( v1->dim != A->m || v2->dim != A-> n )
	error(E_SIZES,"zmv_mltadd");
    
    tracecatch(out = zv_copy(v1,out),"zmv_mltadd");
    
    v2_ve = v2->ve;	out_ve = out->ve;
    m = A->m;	n = A->n;
    
    if ( alpha.re == 0.0 && alpha.im == 0.0 )
	return out;
    
    for ( i = 0; i < m; i++ )
    {
	tmp = __zip__(A->me[i],v2_ve,(int)n,Z_NOCONJ);
	out_ve[i].re += alpha.re*tmp.re - alpha.im*tmp.im;
	out_ve[i].im += alpha.re*tmp.im + alpha.im*tmp.re;
	/**************************************************
	  A_e = A->me[i];
	  sum = 0.0;
	  for ( j = 0; j < n; j++ )
	  sum += A_e[j]*v2_ve[j];
	  out_ve[i] = v1->ve[i] + alpha*sum;
	  **************************************************/
    }
    
    return out;
}

/* zvm_mltadd -- vector-matrix multiply and add a la zvm_mlt()
   -- may not be in situ
   -- returns out == v1 + v2*.A */
ZVEC	*zvm_mltadd(v1,v2,A,alpha,out)
ZVEC	*v1, *v2, *out;
ZMAT	*A;
complex	alpha;
{
    int	/* i, */ j, m, n;
    complex	tmp, /* *A_e, */ *out_ve;
    
    if ( ! v1 || ! v2 || ! A )
	error(E_NULL,"zvm_mltadd");
    if ( v2 == out )
	error(E_INSITU,"zvm_mltadd");
    if ( v1->dim != A->n || A->m != v2->dim )
	error(E_SIZES,"zvm_mltadd");
    
    tracecatch(out = zv_copy(v1,out),"zvm_mltadd");
    
    out_ve = out->ve;	m = A->m;	n = A->n;
    for ( j = 0; j < m; j++ )
    {
	/* tmp = zmlt(v2->ve[j],alpha); */
	tmp.re =   v2->ve[j].re*alpha.re - v2->ve[j].im*alpha.im;
	tmp.im =   v2->ve[j].re*alpha.im + v2->ve[j].im*alpha.re;
	if ( tmp.re != 0.0 || tmp.im != 0.0 )
	    __zmltadd__(out_ve,A->me[j],tmp,(int)n,Z_CONJ);
	/**************************************************
	  A_e = A->me[j];
	  for ( i = 0; i < n; i++ )
	  out_ve[i] += A_e[i]*tmp;
	**************************************************/
    }
    
    return out;
}

/* zget_col -- gets a specified column of a matrix; returned as a vector */
ZVEC	*zget_col(mat,col,vec)
int	col;
ZMAT	*mat;
ZVEC	*vec;
{
	unsigned int	i;

	if ( mat==ZMNULL )
		error(E_NULL,"zget_col");
	if ( col < 0 || col >= mat->n )
		error(E_RANGE,"zget_col");
	if ( vec==ZVNULL || vec->dim<mat->m )
		vec = zv_resize(vec,mat->m);

	for ( i=0; i<mat->m; i++ )
	    vec->ve[i] = mat->me[i][col];

	return (vec);
}

/* zget_row -- gets a specified row of a matrix and retruns it as a vector */
ZVEC	*zget_row(mat,row,vec)
int	row;
ZMAT	*mat;
ZVEC	*vec;
{
	int	/* i, */ lim;

	if ( mat==ZMNULL )
		error(E_NULL,"zget_row");
	if ( row < 0 || row >= mat->m )
		error(E_RANGE,"zget_row");
	if ( vec==ZVNULL || vec->dim<mat->n )
		vec = zv_resize(vec,mat->n);

	lim = min(mat->n,vec->dim);

	/* for ( i=0; i<mat->n; i++ ) */
	/*     vec->ve[i] = mat->me[row][i]; */
	MEMCOPY(mat->me[row],vec->ve,lim,complex);

	return (vec);
}

/* zset_col -- sets column of matrix to values given in vec (in situ) */
ZMAT	*zset_col(mat,col,vec)
ZMAT	*mat;
ZVEC	*vec;
int	col;
{
	unsigned int	i,lim;

	if ( mat==ZMNULL || vec==ZVNULL )
		error(E_NULL,"zset_col");
	if ( col < 0 || col >= mat->n )
		error(E_RANGE,"zset_col");
	lim = min(mat->m,vec->dim);
	for ( i=0; i<lim; i++ )
	    mat->me[i][col] = vec->ve[i];

	return (mat);
}

/* zset_row -- sets row of matrix to values given in vec (in situ) */
ZMAT	*zset_row(mat,row,vec)
ZMAT	*mat;
ZVEC	*vec;
int	row;
{
	unsigned int	/* j, */ lim;

	if ( mat==ZMNULL || vec==ZVNULL )
		error(E_NULL,"zset_row");
	if ( row < 0 || row >= mat->m )
		error(E_RANGE,"zset_row");
	lim = min(mat->n,vec->dim);
	/* for ( j=j0; j<lim; j++ ) */
	/*     mat->me[row][j] = vec->ve[j]; */
	MEMCOPY(vec->ve,mat->me[row],lim,complex);

	return (mat);
}

/* zm_rand -- randomise a complex matrix; uniform in [0,1)+[0,1)*i */
ZMAT	*zm_rand(A)
ZMAT	*A;
{
    int		i;

    if ( ! A )
	error(E_NULL,"zm_rand");

    for ( i = 0; i < A->m; i++ )
	mrandlist((Real *)(A->me[i]),2*A->n);

    return A;
}