zmatop.c

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

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**			     Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/



#include	<stdio.h>
#include	"zmatrix.h"

static	char	rcsid[] = "$Id: zmatop.c,v 1.2 1995/03/27 15:49:03 des Exp $";


#define	is_zero(z)	((z).re == 0.0 && (z).im == 0.0)

/* zm_add -- matrix addition -- may be in-situ */
ZMAT	*zm_add(mat1,mat2,out)
ZMAT	*mat1,*mat2,*out;
{
    unsigned int	m,n,i;
    
    if ( mat1==ZMNULL || mat2==ZMNULL )
	error(E_NULL,"zm_add");
    if ( mat1->m != mat2->m || mat1->n != mat2->n )
	error(E_SIZES,"zm_add");
    if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
	out = zm_resize(out,mat1->m,mat1->n);
    m = mat1->m;	n = mat1->n;
    for ( i=0; i<m; i++ )
    {
	__zadd__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
	/**************************************************
	  for ( j=0; j<n; j++ )
	  out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
	  **************************************************/
    }
    
    return (out);
}

/* zm_sub -- matrix subtraction -- may be in-situ */
ZMAT	*zm_sub(mat1,mat2,out)
ZMAT	*mat1,*mat2,*out;
{
    unsigned int	m,n,i;
    
    if ( mat1==ZMNULL || mat2==ZMNULL )
	error(E_NULL,"zm_sub");
    if ( mat1->m != mat2->m || mat1->n != mat2->n )
	error(E_SIZES,"zm_sub");
    if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
	out = zm_resize(out,mat1->m,mat1->n);
    m = mat1->m;	n = mat1->n;
    for ( i=0; i<m; i++ )
    {
	__zsub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
	/**************************************************
	  for ( j=0; j<n; j++ )
	  out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
	**************************************************/
    }
    
    return (out);
}

/*
  Note: In the following routines, "adjoint" means complex conjugate
  transpose:
  A* = conjugate(A^T)
  */

/* zm_mlt -- matrix-matrix multiplication */
ZMAT	*zm_mlt(A,B,OUT)
ZMAT	*A,*B,*OUT;
{
    unsigned int	i, /* j, */ k, m, n, p;
    complex	**A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
    
    if ( A==ZMNULL || B==ZMNULL )
	error(E_NULL,"zm_mlt");
    if ( A->n != B->m )
	error(E_SIZES,"zm_mlt");
    if ( A == OUT || B == OUT )
	error(E_INSITU,"zm_mlt");
    m = A->m;	n = A->n;	p = B->n;
    A_v = A->me;		B_v = B->me;
    
    if ( OUT==ZMNULL || OUT->m != A->m || OUT->n != B->n )
	OUT = zm_resize(OUT,A->m,B->n);
    
    /****************************************************************
      for ( i=0; i<m; i++ )
      for  ( j=0; j<p; j++ )
      {
      sum = 0.0;
      for ( k=0; k<n; k++ )
      sum += A_v[i][k]*B_v[k][j];
      OUT->me[i][j] = sum;
      }
    ****************************************************************/
    zm_zero(OUT);
    for ( i=0; i<m; i++ )
	for ( k=0; k<n; k++ )
	{
	    if ( ! is_zero(A_v[i][k]) )
		__zmltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p,Z_NOCONJ);
	    /**************************************************
	      B_row = B_v[k];	OUT_row = OUT->me[i];
	      for ( j=0; j<p; j++ )
	      (*OUT_row++) += tmp*(*B_row++);
	    **************************************************/
	}
    
    return OUT;
}

/* zmma_mlt -- matrix-matrix adjoint multiplication
   -- A.B* is returned, and stored in OUT */
ZMAT	*zmma_mlt(A,B,OUT)
ZMAT	*A, *B, *OUT;
{
    int	i, j, limit;
    /* complex	*A_row, *B_row, sum; */
    
    if ( ! A || ! B )
	error(E_NULL,"zmma_mlt");
    if ( A == OUT || B == OUT )
	error(E_INSITU,"zmma_mlt");
    if ( A->n != B->n )
	error(E_SIZES,"zmma_mlt");
    if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
	OUT = zm_resize(OUT,A->m,B->m);
    
    limit = A->n;
    for ( i = 0; i < A->m; i++ )
	for ( j = 0; j < B->m; j++ )
	{
	    OUT->me[i][j] = __zip__(B->me[j],A->me[i],(int)limit,Z_CONJ);
	    /**************************************************
	      sum = 0.0;
	      A_row = A->me[i];
	      B_row = B->me[j];
	      for ( k = 0; k < limit; k++ )
	      sum += (*A_row++)*(*B_row++);
	      OUT->me[i][j] = sum;
	      **************************************************/
	}
    
    return OUT;
}

/* zmam_mlt -- matrix adjoint-matrix multiplication
   -- A*.B is returned, result stored in OUT */
ZMAT	*zmam_mlt(A,B,OUT)
ZMAT	*A, *B, *OUT;
{
    int	i, k, limit;
    /* complex	*B_row, *OUT_row, multiplier; */
    complex	tmp;
    
    if ( ! A || ! B )
	error(E_NULL,"zmam_mlt");
    if ( A == OUT || B == OUT )
	error(E_INSITU,"zmam_mlt");
    if ( A->m != B->m )
	error(E_SIZES,"zmam_mlt");
    if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
	OUT = zm_resize(OUT,A->n,B->n);
    
    limit = B->n;
    zm_zero(OUT);
    for ( k = 0; k < A->m; k++ )
	for ( i = 0; i < A->n; i++ )
	{
	    tmp.re =   A->me[k][i].re;
	    tmp.im = - A->me[k][i].im;
	    if ( ! is_zero(tmp) )
		__zmltadd__(OUT->me[i],B->me[k],tmp,(int)limit,Z_NOCONJ);
	}
    
    return OUT;
}

/* zmv_mlt -- matrix-vector multiplication 
   -- Note: b is treated as a column vector */
ZVEC	*zmv_mlt(A,b,out)
ZMAT	*A;
ZVEC	*b,*out;
{
    unsigned int	i, m, n;
    complex	**A_v, *b_v /*, *A_row */;
    /* register complex	sum; */
    
    if ( A==ZMNULL || b==ZVNULL )
	error(E_NULL,"zmv_mlt");
    if ( A->n != b->dim )
	error(E_SIZES,"zmv_mlt");
    if ( b == out )
	error(E_INSITU,"zmv_mlt");
    if ( out == ZVNULL || out->dim != A->m )
	out = zv_resize(out,A->m);
    
    m = A->m;		n = A->n;
    A_v = A->me;		b_v = b->ve;
    for ( i=0; i<m; i++ )
    {
	/* for ( j=0; j<n; j++ )
	   sum += A_v[i][j]*b_v[j]; */
	out->ve[i] = __zip__(A_v[i],b_v,(int)n,Z_NOCONJ);
	/**************************************************
	  A_row = A_v[i];		b_v = b->ve;
	  for ( j=0; j<n; j++ )
	  sum += (*A_row++)*(*b_v++);
	  out->ve[i] = sum;
	**************************************************/
    }
    
    return out;
}

/* zsm_mlt -- scalar-matrix multiply -- may be in-situ */
ZMAT	*zsm_mlt(scalar,matrix,out)
complex	scalar;
ZMAT	*matrix,*out;
{
    unsigned int	m,n,i;
    
    if ( matrix==ZMNULL )
	error(E_NULL,"zsm_mlt");
    if ( out==ZMNULL || out->m != matrix->m || out->n != matrix->n )
	out = zm_resize(out,matrix->m,matrix->n);
    m = matrix->m;	n = matrix->n;
    for ( i=0; i<m; i++ )
	__zmlt__(matrix->me[i],scalar,out->me[i],(int)n);
    /**************************************************
      for ( j=0; j<n; j++ )
      out->me[i][j] = scalar*matrix->me[i][j];
      **************************************************/
    return (out);
}

/* zvm_mlt -- vector adjoint-matrix multiplication */
ZVEC	*zvm_mlt(A,b,out)
ZMAT	*A;
ZVEC	*b,*out;
{
    unsigned int	j,m,n;
    /* complex	sum,**A_v,*b_v; */
    
    if ( A==ZMNULL || b==ZVNULL )
	error(E_NULL,"zvm_mlt");
    if ( A->m != b->dim )
	error(E_SIZES,"zvm_mlt");
    if ( b == out )
	error(E_INSITU,"zvm_mlt");
    if ( out == ZVNULL || out->dim != A->n )
	out = zv_resize(out,A->n);
    
    m = A->m;		n = A->n;
    
    zv_zero(out);
    for ( j = 0; j < m; j++ )
	if ( b->ve[j].re != 0.0 || b->ve[j].im != 0.0  )
	    __zmltadd__(out->ve,A->me[j],b->ve[j],(int)n,Z_CONJ);
    /**************************************************
      A_v = A->me;		b_v = b->ve;
      for ( j=0; j<n; j++ )
      {
      sum = 0.0;
      for ( i=0; i<m; i++ )
      sum += b_v[i]*A_v[i][j];
      out->ve[j] = sum;
      }
      **************************************************/
    
    return out;
}

/* zm_adjoint -- adjoint matrix */
ZMAT	*zm_adjoint(in,out)
ZMAT	*in, *out;
{
    int	i, j;
    int	in_situ;
    complex	tmp;
    
    if ( in == ZMNULL )
	error(E_NULL,"zm_adjoint");