bdfactor.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.
**
***************************************************************************/


/*
  Band matrix factorisation routines
  */

/* bdfactor.c  18/11/93 */
static	char	rcsid[] = "$Id: ";

#include	<stdio.h>
#include	<math.h>
#include        "matrix2.h"


/* generate band matrix 
   for a matrix  with n columns,
   lb subdiagonals and ub superdiagonals;

   Way of saving a band of a matrix:
   first we save subdiagonals (from 0 to lb-1);
   then main diagonal (in the lb row)
   and then superdiagonals (from lb+1 to lb+ub)
   in such a way that the elements which were previously
   in one column are now also in one column
*/
#ifndef ANSI_C
BAND *bd_get(lb,ub,n)
int lb, ub, n;
#else
BAND *bd_get(int lb, int ub, int n)
#endif
{
   BAND *A;

   if (lb < 0 || ub < 0 || n <= 0)
     error(E_NEG,"bd_get");

   if ((A = NEW(BAND)) == (BAND *)NULL)
     error(E_MEM,"bd_get");
   else if (mem_info_is_on()) {
      mem_bytes(TYPE_BAND,0,sizeof(BAND));
      mem_numvar(TYPE_BAND,1);
   }

   lb = A->lb = min(n-1,lb);
   ub = A->ub = min(n-1,ub);
   A->mat = m_get(lb+ub+1,n);
   return A;
}

/* bd_free -- frees BAND matrix -- returns (-1) on error and 0 otherwise */
#ifndef ANSI_C
int bd_free(A)
BAND *A;
#else
int bd_free(BAND *A)
#endif
{
   if ( A == (BAND *)NULL || A->lb < 0 || A->ub < 0 )
     /* don't trust it */
     return (-1);

   if (A->mat) m_free(A->mat);

   if (mem_info_is_on()) {
      mem_bytes(TYPE_BAND,sizeof(BAND),0);
      mem_numvar(TYPE_BAND,-1);
   }

   free((char *)A);
   return 0;
}


/* resize band matrix */
#ifndef ANSI_C
BAND *bd_resize(A,new_lb,new_ub,new_n)
BAND *A;
int new_lb,new_ub,new_n;
#else
BAND *bd_resize(BAND *A, int new_lb, int new_ub, int new_n)
#endif
{
   int lb,ub,i,j,l,shift,umin;
   Real **Av;

   if (new_lb < 0 || new_ub < 0 || new_n <= 0)
     error(E_NEG,"bd_resize");
   if ( ! A )
     return bd_get(new_lb,new_ub,new_n);
    if ( A->lb+A->ub+1 > A->mat->m )
	error(E_INTERN,"bd_resize");

   if ( A->lb == new_lb && A->ub == new_ub && A->mat->n == new_n )
	return A;

   lb = A->lb;
   ub = A->ub;
   Av = A->mat->me;
   umin = min(ub,new_ub);

    /* ensure that unused triangles at edges are zero'd */

   for ( i = 0; i < lb; i++ )
      for ( j = A->mat->n - lb + i; j < A->mat->n; j++ )
	Av[i][j] = 0.0;  
    for ( i = lb+1,l=1; l <= umin; i++,l++ )
      for ( j = 0; j < l; j++ )
	Av[i][j] = 0.0; 

   new_lb = A->lb = min(new_lb,new_n-1);
   new_ub = A->ub = min(new_ub,new_n-1);
   A->mat = m_resize(A->mat,new_lb+new_ub+1,new_n);
   Av = A->mat->me;

   /* if new_lb != lb then move the rows to get the main diag 
      in the new_lb row */

   if (new_lb > lb) {
      shift = new_lb-lb;

      for (i=lb+umin, l=i+shift; i >= 0; i--,l--)
	MEM_COPY(Av[i],Av[l],new_n*sizeof(Real));
      for (l=shift-1; l >= 0; l--)
	__zero__(Av[l],new_n);
   }
   else if (new_lb < lb) { 
      shift = lb - new_lb;

      for (i=shift, l=0; i <= lb+umin; i++,l++)
	MEM_COPY(Av[i],Av[l],new_n*sizeof(Real));
      for (i=lb+umin+1; i <= new_lb+new_ub; i++)
	__zero__(Av[i],new_n);
   }

   return A;
}


/* bd_copy -- copies band matrix A to B, returning B
	-- if B is NULL, create
	-- B is set to the correct size */
#ifndef ANSI_C
BAND *bd_copy(A,B)
BAND *A,*B;
#else
BAND *bd_copy(const BAND *A, BAND *B)
#endif
{
   int lb,ub,i,j,n;
   
   if ( !A )
     error(E_NULL,"bd_copy");

   if (A == B) return B;
   
   n = A->mat->n;
   if ( !B )
     B = bd_get(A->lb,A->ub,n);
   else if (B->lb != A->lb || B->ub != A->ub || B->mat->n != n )
     B = bd_resize(B,A->lb,A->ub,n);
   
   if (A->mat == B->mat) return B;
   ub = B->ub = A->ub;
   lb = B->lb = A->lb;

   for ( i=0, j=n-lb; i <= lb; i++, j++ )
     MEM_COPY(A->mat->me[i],B->mat->me[i],j*sizeof(Real));   

   for ( i=lb+1, j=1; i <= lb+ub; i++, j++ )
     MEM_COPY(A->mat->me[i]+j,B->mat->me[i]+j,(n - j)*sizeof(Real));     

   return B;
}


/* copy band matrix bA to a square matrix A returning A */
#ifndef ANSI_C
MAT *band2mat(bA,A)
BAND *bA;
MAT *A;
#else
MAT *band2mat(const BAND *bA, MAT *A)
#endif
{
   int i,j,l,n,n1;
   int lb, ub;
   Real **bmat;

   if ( !bA )
     error(E_NULL,"band2mat");
   if ( bA->mat == A )
     error(E_INSITU,"band2mat");

   ub = bA->ub;
   lb = bA->lb;
   n = bA->mat->n;
   n1 = n-1;
   bmat = bA->mat->me;

   A = m_resize(A,n,n);
   m_zero(A);

   for (j=0; j < n; j++)
     for (i=min(n1,j+lb),l=lb+j-i; i >= max(0,j-ub); i--,l++)
       A->me[i][j] = bmat[l][j];

   return A;
}

/* copy a square matrix to a band matrix with 
   lb subdiagonals and ub superdiagonals */
#ifndef ANSI_C
BAND *mat2band(A,lb,ub,bA)
BAND *bA;
MAT *A;
int lb, ub;
#else
BAND *mat2band(const MAT *A, int lb, int ub,BAND *bA)
#endif
{
   int i, j, l, n1;
   Real **bmat;
   
   if (! A )
     error(E_NULL,"mat2band");
   if (ub < 0 || lb < 0)
     error(E_SIZES,"mat2band");
   if ( bA != (BAND *)NULL && bA->mat == A )
     error(E_INSITU,"mat2band");

   n1 = A->n-1;
   lb = min(n1,lb);
   ub = min(n1,ub);
   bA = bd_resize(bA,lb,ub,n1+1);
   bmat = bA->mat->me;

   for (j=0; j <= n1; j++)
     for (i=min(n1,j+lb),l=lb+j-i; i >= max(0,j-ub); i--,l++)
       bmat[l][j] = A->me[i][j];

   return bA;
}



/* transposition of matrix in;
   out - matrix after transposition;
   can be done in situ
*/
#ifndef ANSI_C
BAND *bd_transp(in,out)
BAND *in, *out;
#else
BAND *bd_transp(const BAND *in, BAND *out)
#endif
{
   int i, j, jj, l, k, lb, ub, lub, n, n1;
   int in_situ;
   Real  **in_v, **out_v;
   
   if ( in == (BAND *)NULL || in->mat == (MAT *)NULL )
     error(E_NULL,"bd_transp");

   lb = in->lb;
   ub = in->ub;
   lub = lb+ub;
   n = in->mat->n;
   n1 = n-1;

   in_situ = ( in == out );
   if ( ! in_situ )
       out = bd_resize(out,ub,lb,n);
   else
   {   /* only need to swap lb and ub fields */
       out->lb = ub;
       out->ub = lb;
   }

   in_v = in->mat->me;
   
   if (! in_situ) {
      int sh_in,sh_out; 

      out_v = out->mat->me;
      for (i=0, l=lub, k=lb-i; i <= lub; i++,l--,k--) {
	 sh_in = max(-k,0);
	 sh_out = max(k,0);
	 MEM_COPY(&(in_v[i][sh_in]),&(out_v[l][sh_out]),
		  (n-sh_in-sh_out)*sizeof(Real));
	 /**********************************
	 for (j=n1-sh_out, jj=n1-sh_in; j >= sh_in; j--,jj--) {
	    out_v[l][jj] = in_v[i][j];
	 }
	 **********************************/
      }
   }
   else if (ub == lb) {
      Real tmp;

      for (i=0, l=lub, k=lb-i; i < lb; i++,l--,k--) {
	 for (j=n1-k, jj=n1; j >= 0; j--,jj--) {
	    tmp = in_v[l][jj];
	    in_v[l][jj] = in_v[i][j];
	    in_v[i][j] = tmp;
	 }
      }
   }
   else if (ub > lb) {  /* hence i-ub <= 0 & l-lb >= 0 */
      int p,pp,lbi;
      
      for (i=0, l=lub; i < (lub+1)/2; i++,l--) {
	 lbi = lb-i;
	 for (j=l-lb, jj=0, p=max(-lbi,0), pp = max(l-ub,0); j <= n1; 
	      j++,jj++,p++,pp++) {
	    in_v[l][pp] = in_v[i][p];
	    in_v[i][jj] = in_v[l][j];
	 }
	 for (  ; p <= n1-max(lbi,0); p++,pp++)
	   in_v[l][pp] = in_v[i][p];
      }
      
      if (lub%2 == 0) { /* shift only */
	 i = lub/2;
	 for (j=max(i-lb,0), jj=0; jj <= n1-ub+i; j++,jj++) 
	   in_v[i][jj] = in_v[i][j];
      }
   }
   else {      /* ub < lb, hence ub-l <= 0 & lb-i >= 0 */
      int p,pp,ubi;

      for (i=0, l=lub; i < (lub+1)/2; i++,l--) {
	 ubi = i-ub;
	 for (j=n1-max(lb-l,0), jj=n1-max(-ubi,0), p=n1-lb+i, pp=n1;
	      p >= 0; j--, jj--, pp--, p--) {
	    in_v[i][jj] = in_v[l][j];
	    in_v[l][pp] = in_v[i][p];
	 }
	 for (  ; jj >= max(ubi,0); j--, jj--)
	   in_v[i][jj] = in_v[l][j];
      }

      if (lub%2 == 0) {  /* shift only */
	 i = lub/2;
	 for (j=n1-lb+i, jj=n1-max(ub-i,0); j >= 0; j--, jj--) 
	    in_v[i][jj] = in_v[i][j];
      }
   }

   return out;
}

/* bdv_mltadd -- band matrix-vector multiply and add
   -- returns out <- x + s.bA.y
   -- if y is NULL then create y (as zero vector)
   -- error if either A or x is NULL */
#ifndef ANSI_C
VEC	*bdv_mltadd(x,y,bA,s,out)
     BAND	*bA;
     VEC	*x, *y;
     double	s;
     VEC *out;
#else
VEC	*bdv_mltadd(const VEC *x, const VEC *y, const BAND *bA,
		    double s, VEC *out)
#endif
{
  int	i, j;

  if ( ! bA || ! x || ! y )
    error(E_NULL,"bdv_mltadd");
  if ( bA->mat->n != x->dim || y->dim != x->dim )
    error(E_SIZES,"bdv_mltadd");
  if ( ! out || out->dim != x->dim )
    out = v_resize(out,x->dim);
  out = v_copy(x,out);

  for ( j = 0; j < x->dim; j++ )
    for ( i = max(j-bA->ub,0); i <= j+bA->lb && i < x->dim; i++ )
      out->ve[i] += s*bd_get_val(bA,i,j)*y->ve[j];

  return out;
}

/* vbd_mltadd -- band matrix-vector multiply and add