iternsym.c

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

/**************************************************************************
**
** Copyright (C) 1993 David E. Stewart & 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.
**
***************************************************************************/


/* iter.c 17/09/93 */

/* 
  ITERATIVE METHODS - implementation of several iterative methods;
  see also iter0.c
*/

#include        <stdio.h>
#include	<math.h>
#include        "matrix.h"
#include        "matrix2.h"
#include	"sparse.h"
#include        "iter.h"

static char rcsid[] = "$Header: iternsym.c,v 1.6 1995/01/30 14:53:01 des Exp $";


#ifdef ANSI_C
VEC	*spCHsolve(SPMAT *,VEC *,VEC *);
#else
VEC	*spCHsolve();
#endif


/* 
  iter_cgs -- uses CGS to compute a solution x to A.x=b
*/
#ifndef ANSI_C
VEC	*iter_cgs(ip,r0)
ITER *ip;
VEC *r0;
#else
VEC	*iter_cgs(ITER *ip, VEC *r0)
#endif
{
   STATIC VEC  *p = VNULL, *q = VNULL, *r = VNULL, *u = VNULL;
   STATIC VEC  *v = VNULL, *z = VNULL;
   VEC  *tmp;
   Real	alpha, beta, nres, rho, old_rho, sigma, inner;

   if (ip == INULL)
     error(E_NULL,"iter_cgs");
   if (!ip->Ax || !ip->b || !r0)
     error(E_NULL,"iter_cgs");
   if ( ip->x == ip->b )
     error(E_INSITU,"iter_cgs");
   if (!ip->stop_crit)
     error(E_NULL,"iter_cgs");
   if ( r0->dim != ip->b->dim )
     error(E_SIZES,"iter_cgs");
   
   if ( ip->eps <= 0.0 ) ip->eps = MACHEPS;
   
   p = v_resize(p,ip->b->dim);
   q = v_resize(q,ip->b->dim);
   r = v_resize(r,ip->b->dim);
   u = v_resize(u,ip->b->dim);
   v = v_resize(v,ip->b->dim);

   MEM_STAT_REG(p,TYPE_VEC);
   MEM_STAT_REG(q,TYPE_VEC);
   MEM_STAT_REG(r,TYPE_VEC);
   MEM_STAT_REG(u,TYPE_VEC);
   MEM_STAT_REG(v,TYPE_VEC);

   if (ip->Bx) {
      z = v_resize(z,ip->b->dim);
      MEM_STAT_REG(z,TYPE_VEC); 
   }

   if (ip->x != VNULL) {
      if (ip->x->dim != ip->b->dim)
	error(E_SIZES,"iter_cgs");
      ip->Ax(ip->A_par,ip->x,v);    		/* v = A*x */
      if (ip->Bx) {
	 v_sub(ip->b,v,v);			/* v = b - A*x */
	 (ip->Bx)(ip->B_par,v,r);		/* r = B*(b-A*x) */
      }
      else v_sub(ip->b,v,r);			/* r = b-A*x */
   }
   else {  /* ip->x == 0 */
      ip->x = v_get(ip->b->dim);		/* x == 0 */
      ip->shared_x = FALSE;
      if (ip->Bx) (ip->Bx)(ip->B_par,ip->b,r);    /* r = B*b */
      else v_copy(ip->b,r);                       /* r = b */
   }

   v_zero(p);	
   v_zero(q);
   old_rho = 1.0;
   
   for (ip->steps = 0; ip->steps <= ip->limit; ip->steps++) {

      inner = in_prod(r,r);
      nres = sqrt(fabs(inner));
      if (ip->steps == 0) ip->init_res = nres;

      if (ip->info) ip->info(ip,nres,r,VNULL);
      if ( ip->stop_crit(ip,nres,r,VNULL) ) break;

      rho = in_prod(r0,r);
      if ( old_rho == 0.0 )
	error(E_BREAKDOWN,"iter_cgs");
      beta = rho/old_rho;
      v_mltadd(r,q,beta,u);
      v_mltadd(q,p,beta,v);
      v_mltadd(u,v,beta,p);
      
      (ip->Ax)(ip->A_par,p,q);
      if (ip->Bx) {
	 (ip->Bx)(ip->B_par,q,z);
	 tmp = z;
      }
      else tmp = q;
      
      sigma = in_prod(r0,tmp);
      if ( sigma == 0.0 )
	error(E_BREAKDOWN,"iter_cgs");
      alpha = rho/sigma;
      v_mltadd(u,tmp,-alpha,q);
      v_add(u,q,v);
      
      (ip->Ax)(ip->A_par,v,u);
      if (ip->Bx) {
	 (ip->Bx)(ip->B_par,u,z);
	 tmp = z;
      }
      else tmp = u;
      
      v_mltadd(r,tmp,-alpha,r);
      v_mltadd(ip->x,v,alpha,ip->x);
      
      old_rho = rho;
   }

#ifdef THREADSAFE
   V_FREE(p);	V_FREE(q);	V_FREE(r);	V_FREE(u);
   V_FREE(v);	V_FREE(z);
#endif

   return ip->x;
}



/* iter_spcgs -- simple interface for SPMAT data structures 
   use always as follows:
      x = iter_spcgs(A,B,b,r0,tol,x,limit,steps);
   or 
      x = iter_spcgs(A,B,b,r0,tol,VNULL,limit,steps);
   In the second case the solution vector is created.  
   If B is not NULL then it is a preconditioner. 
*/
#ifndef ANSI_C
VEC	*iter_spcgs(A,B,b,r0,tol,x,limit,steps)
SPMAT	*A, *B;
VEC	*b, *r0, *x;
double	tol;
int     *steps,limit;
#else
VEC	*iter_spcgs(SPMAT *A, SPMAT *B, VEC *b, VEC *r0, double tol,
		    VEC *x, int limit, int *steps)
#endif
{	
   ITER *ip;
   
   ip = iter_get(0,0);
   ip->Ax = (Fun_Ax) sp_mv_mlt;
   ip->A_par = (void *) A;
   if (B) {
      ip->Bx = (Fun_Ax) sp_mv_mlt;
      ip->B_par = (void *) B;
   }
   else {
      ip->Bx = (Fun_Ax) NULL;
      ip->B_par = NULL;
   }
   ip->info = (Fun_info) NULL;
   ip->limit = limit;
   ip->b = b;
   ip->eps = tol;
   ip->x = x;
   iter_cgs(ip,r0);
   x = ip->x;
   if (steps) *steps = ip->steps;
   ip->shared_x = ip->shared_b = TRUE;   
   iter_free(ip);   /* release only ITER structure */
   return x;		

}

/*
  Routine for performing LSQR -- the least squares QR algorithm
  of Paige and Saunders:
  "LSQR: an algorithm for sparse linear equations and
  sparse least squares", ACM Trans. Math. Soft., v. 8
  pp. 43--71 (1982)
  */
/* iter_lsqr -- sparse CG-like least squares routine:
   -- finds min_x ||A.x-b||_2 using A defined through A & AT
   -- returns x (if x != NULL) */
#ifndef ANSI_C
VEC	*iter_lsqr(ip)
ITER *ip;
#else
VEC	*iter_lsqr(ITER *ip)
#endif
{
   STATIC VEC	*u = VNULL, *v = VNULL, *w = VNULL, *tmp = VNULL;
   Real	alpha, beta, phi, phi_bar;
   Real rho, rho_bar, rho_max, theta, nres;
   Real	s, c;	/* for Givens' rotations */
   int  m, n;
   
   if ( ! ip || ! ip->b || !ip->Ax || !ip->ATx )
     error(E_NULL,"iter_lsqr");
   if ( ip->x == ip->b )
     error(E_INSITU,"iter_lsqr");
   if (!ip->stop_crit || !ip->x)
     error(E_NULL,"iter_lsqr");

   if ( ip->eps <= 0.0 ) ip->eps = MACHEPS;
   
   m = ip->b->dim;	
   n = ip->x->dim;

   u = v_resize(u,(unsigned int)m);
   v = v_resize(v,(unsigned int)n);
   w = v_resize(w,(unsigned int)n);
   tmp = v_resize(tmp,(unsigned int)n);

   MEM_STAT_REG(u,TYPE_VEC);
   MEM_STAT_REG(v,TYPE_VEC);
   MEM_STAT_REG(w,TYPE_VEC);
   MEM_STAT_REG(tmp,TYPE_VEC);  

   if (ip->x != VNULL) {
      ip->Ax(ip->A_par,ip->x,u);    		/* u = A*x */
      v_sub(ip->b,u,u);				/* u = b-A*x */
   }
   else {  /* ip->x == 0 */
      ip->x = v_get(ip->b->dim);
      ip->shared_x = FALSE;
      v_copy(ip->b,u);                       /* u = b */
   }
 
   beta = v_norm2(u); 
   if ( beta == 0.0 ) return ip->x;

   sv_mlt(1.0/beta,u,u);
   (ip->ATx)(ip->AT_par,u,v);
   alpha = v_norm2(v);
   if ( alpha == 0.0 ) return ip->x;

   sv_mlt(1.0/alpha,v,v);
   v_copy(v,w);
   phi_bar = beta;
   rho_bar = alpha;
   
   rho_max = 1.0;
   for (ip->steps = 0; ip->steps <= ip->limit; ip->steps++) {

      tmp = v_resize(tmp,m);
      (ip->Ax)(ip->A_par,v,tmp);
      
      v_mltadd(tmp,u,-alpha,u);
      beta = v_norm2(u);	
      sv_mlt(1.0/beta,u,u);
      
      tmp = v_resize(tmp,n);
      (ip->ATx)(ip->AT_par,u,tmp);
      v_mltadd(tmp,v,-beta,v);
      alpha = v_norm2(v);	
      sv_mlt(1.0/alpha,v,v);
      
      rho = sqrt(rho_bar*rho_bar+beta*beta);
      if ( rho > rho_max )
	rho_max = rho;
      c   = rho_bar/rho;
      s   = beta/rho;
      theta   =  s*alpha;
      rho_bar = -c*alpha;
      phi     =  c*phi_bar;
      phi_bar =  s*phi_bar;
      
      /* update ip->x & w */
      if ( rho == 0.0 )
	error(E_BREAKDOWN,"iter_lsqr");
      v_mltadd(ip->x,w,phi/rho,ip->x);
      v_mltadd(v,w,-theta/rho,w);

      nres = fabs(phi_bar*alpha*c)*rho_max;

      if (ip->info) ip->info(ip,nres,w,VNULL);
      if (ip->steps == 0) ip->init_res = nres;
      if ( ip->stop_crit(ip,nres,w,VNULL) ) break;
   } 

#ifdef THREADSAFE
   V_FREE(u);	V_FREE(v);	V_FREE(w);	V_FREE(tmp);
#endif

   return ip->x;
}

/* iter_splsqr -- simple interface for SPMAT data structures */
#ifndef ANSI_C
VEC	*iter_splsqr(A,b,tol,x,limit,steps)
SPMAT	*A;
VEC	*b, *x;
double	tol;
int *steps,limit;
#else
VEC	*iter_splsqr(SPMAT *A, VEC *b, double tol, 
		     VEC *x, int limit, int *steps)
#endif
{
   ITER *ip;
   
   ip = iter_get(0,0);
   ip->Ax = (Fun_Ax) sp_mv_mlt;
   ip->A_par = (void *) A;
   ip->ATx = (Fun_Ax) sp_vm_mlt;
   ip->AT_par = (void *) A;
   ip->Bx = (Fun_Ax) NULL;
   ip->B_par = NULL;

   ip->info = (Fun_info) NULL;
   ip->limit = limit;
   ip->b = b;
   ip->eps = tol;
   ip->x = x;
   iter_lsqr(ip);
   x = ip->x;
   if (steps) *steps = ip->steps;
   ip->shared_x = ip->shared_b = TRUE;
   iter_free(ip);   /* release only ITER structure */
   return x;		
}



/* iter_arnoldi -- an implementation of the Arnoldi method;
   iterative refinement is applied.
*/
#ifndef ANSI_C
MAT	*iter_arnoldi_iref(ip,h_rem,Q,H)
ITER  *ip;
Real  *h_rem;
MAT   *Q, *H;
#else
MAT	*iter_arnoldi_iref(ITER *ip, Real *h_rem, MAT *Q, MAT *H)
#endif
{
   STATIC VEC *u=VNULL, *r=VNULL, *s=VNULL, *tmp=VNULL;
   VEC v;     /* auxiliary vector */
   int	i,j;
   Real	h_val, c;
   
   if (ip == INULL)
     error(E_NULL,"iter_arnoldi_iref");
   if ( ! ip->Ax || ! Q || ! ip->x )
     error(E_NULL,"iter_arnoldi_iref");
   if ( ip->k <= 0 )
     error(E_BOUNDS,"iter_arnoldi_iref");
   if ( Q->n != ip->x->dim ||	Q->m != ip->k )
     error(E_SIZES,"iter_arnoldi_iref");
   
   m_zero(Q);
   H = m_resize(H,ip->k,ip->k);
   m_zero(H);

   u = v_resize(u,ip->x->dim);
   r = v_resize(r,ip->k);
   s = v_resize(s,ip->k);
   tmp = v_resize(tmp,ip->x->dim);
   MEM_STAT_REG(u,TYPE_VEC);
   MEM_STAT_REG(r,TYPE_VEC);
   MEM_STAT_REG(s,TYPE_VEC);
   MEM_STAT_REG(tmp,TYPE_VEC);

   v.dim = v.max_dim = ip->x->dim;

   c = v_norm2(ip->x);
   if ( c <= 0.0)
     return H;
   else {
      v.ve = Q->me[0];
      sv_mlt(1.0/c,ip->x,&v);
   }

   v_zero(r);
   v_zero(s);
   for ( i = 0; i < ip->k; i++ )
   {
      v.ve = Q->me[i];
      u = (ip->Ax)(ip->A_par,&v,u);
      for (j = 0; j <= i; j++) {
	 v.ve = Q->me[j];
	 /* modified Gram-Schmidt */
	 r->ve[j] = in_prod(&v,u);
	 v_mltadd(u,&v,-r->ve[j],u);
      }
      h_val = v_norm2(u);
      /* if u == 0 then we have an exact subspace */
      if ( h_val <= 0.0 )
      {
	 *h_rem = h_val;
	 return H;
      }
      /* iterative refinement -- ensures near orthogonality */
      do {
	 v_zero(tmp);
	 for (j = 0; j <= i; j++) {
	    v.ve = Q->me[j];
	    s->ve[j] = in_prod(&v,u);
	    v_mltadd(tmp,&v,s->ve[j],tmp);
	 }
	 v_sub(u,tmp,u);
         v_add(r,s,r);
      } while ( v_norm2(s) > 0.1*(h_val = v_norm2(u)) );
      /* now that u is nearly orthogonal to Q, update H */
      set_col(H,i,r);
      /* check once again if h_val is zero */
      if ( h_val <= 0.0 )
      {
	 *h_rem = h_val;
	 return H;
      }
      if ( i == ip->k-1 )
      {
	 *h_rem = h_val;
	 continue;
      }
      /* H->me[i+1][i] = h_val; */
      m_set_val(H,i+1,i,h_val);
      v.ve = Q->me[i+1];
      sv_mlt(1.0/h_val,u,&v);
   }