iternsym.c

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

#ifdef THREADSAFE
   V_FREE(u);   V_FREE(r);   V_FREE(s);   V_FREE(tmp);
#endif

   return H;
}

/* iter_arnoldi -- an implementation of the Arnoldi method;
   modified Gram-Schmidt algorithm
*/
#ifndef ANSI_C
MAT	*iter_arnoldi(ip,h_rem,Q,H)
ITER  *ip;
Real  *h_rem;
MAT   *Q, *H;
#else
MAT	*iter_arnoldi(ITER *ip, Real *h_rem, MAT *Q, MAT *H)
#endif
{
   STATIC VEC *u=VNULL, *r=VNULL;
   VEC v;     /* auxiliary vector */
   int	i,j;
   Real	h_val, c;
   
   if (ip == INULL)
     error(E_NULL,"iter_arnoldi");
   if ( ! ip->Ax || ! Q || ! ip->x )
     error(E_NULL,"iter_arnoldi");
   if ( ip->k <= 0 )
     error(E_BOUNDS,"iter_arnoldi");
   if ( Q->n != ip->x->dim ||	Q->m != ip->k )
     error(E_SIZES,"iter_arnoldi");
   
   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);
   MEM_STAT_REG(u,TYPE_VEC);
   MEM_STAT_REG(r,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);
   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;
      }
      set_col(H,i,r);
      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);
   }

#ifdef THREADSAFE
   V_FREE(u);	V_FREE(r);
#endif
   
   return H;
}



/* iter_sparnoldi -- uses arnoldi() with an explicit representation of A */
#ifndef ANSI_C
MAT	*iter_sparnoldi(A,x0,m,h_rem,Q,H)
SPMAT	*A;
VEC	*x0;
int	m;
Real	*h_rem;
MAT	*Q, *H;
#else
MAT	*iter_sparnoldi(SPMAT *A, VEC *x0, int m, Real *h_rem, MAT *Q, MAT *H)
#endif
{
   ITER *ip;
   
   ip = iter_get(0,0);
   ip->Ax = (Fun_Ax) sp_mv_mlt;
   ip->A_par = (void *) A;
   ip->x = x0;
   ip->k = m;
   iter_arnoldi_iref(ip,h_rem,Q,H);
   ip->shared_x = ip->shared_b = TRUE;
   iter_free(ip);   /* release only ITER structure */
   return H;	
}


/* for testing gmres */
#ifndef ANSI_C
static void test_gmres(ip,i,Q,R,givc,givs,h_val)
ITER *ip;
int i;
MAT *Q, *R;
VEC *givc, *givs;
double h_val;
#else
static void test_gmres(ITER *ip, int i, MAT *Q, MAT *R,
		       VEC *givc, VEC *givs, double h_val)
#endif
{
   VEC vt, vt1;
   STATIC MAT *Q1=MNULL, *R1=MNULL;
   int j;
   
   /* test Q*A*Q^T = R  */

   Q = m_resize(Q,i+1,ip->b->dim);
   Q1 = m_resize(Q1,i+1,ip->b->dim);
   R1 = m_resize(R1,i+1,i+1);
   MEM_STAT_REG(Q1,TYPE_MAT);
   MEM_STAT_REG(R1,TYPE_MAT);

   vt.dim = vt.max_dim = ip->b->dim;
   vt1.dim = vt1.max_dim = ip->b->dim;
   for (j=0; j <= i; j++) {
      vt.ve = Q->me[j];
      vt1.ve = Q1->me[j];
      ip->Ax(ip->A_par,&vt,&vt1);
   }

   mmtr_mlt(Q,Q1,R1);
   R1 = m_resize(R1,i+2,i+1);
   for (j=0; j < i; j++)
     R1->me[i+1][j] = 0.0;
   R1->me[i+1][i] = h_val;
   
   for (j = 0; j <= i; j++) {
      rot_rows(R1,j,j+1,givc->ve[j],givs->ve[j],R1);
   }

   R1 = m_resize(R1,i+1,i+1);
   m_sub(R,R1,R1);
   /* if (m_norm_inf(R1) > MACHEPS*ip->b->dim)  */
#ifndef MEX
   printf(" %d. ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
	  ip->steps,m_norm_inf(R1),MACHEPS);
#endif
   
   /* check Q*Q^T = I */
   
   Q = m_resize(Q,i+1,ip->b->dim);
   mmtr_mlt(Q,Q,R1);
   for (j=0; j <= i; j++)
     R1->me[j][j] -= 1.0;
#ifndef MEX
   if (m_norm_inf(R1) > MACHEPS*ip->b->dim)
     printf(" ! m_norm_inf(Q*Q^T) = %g\n",m_norm_inf(R1));  
#endif
#ifdef THREADSAFE
   M_FREE(Q1);	M_FREE(R1);
#endif
}


/* gmres -- generalised minimum residual algorithm of Saad & Schultz
   SIAM J. Sci. Stat. Comp. v.7, pp.856--869 (1986)
*/
#ifndef ANSI_C
VEC	*iter_gmres(ip)
ITER *ip;
#else
VEC	*iter_gmres(ITER *ip)
#endif
{
   STATIC VEC *u=VNULL, *r=VNULL, *rhs = VNULL;
   STATIC VEC *givs=VNULL, *givc=VNULL, *z = VNULL;
   STATIC MAT *Q = MNULL, *R = MNULL;
   VEC *rr, v, v1;   /* additional pointers (not real vectors) */
   int	i,j, done;
   Real	nres;
/*   Real last_h;  */
   
   if (ip == INULL)
     error(E_NULL,"iter_gmres");
   if ( ! ip->Ax || ! ip->b )
     error(E_NULL,"iter_gmres");
   if ( ! ip->stop_crit )
     error(E_NULL,"iter_gmres");
   if ( ip->k <= 0 )
     error(E_BOUNDS,"iter_gmres");
   if (ip->x != VNULL && ip->x->dim != ip->b->dim)
     error(E_SIZES,"iter_gmres");
   if (ip->eps <= 0.0) ip->eps = MACHEPS;

   r = v_resize(r,ip->k+1);
   u = v_resize(u,ip->b->dim);
   rhs = v_resize(rhs,ip->k+1);
   givs = v_resize(givs,ip->k);  /* Givens rotations */
   givc = v_resize(givc,ip->k); 
   
   MEM_STAT_REG(r,TYPE_VEC);
   MEM_STAT_REG(u,TYPE_VEC);
   MEM_STAT_REG(rhs,TYPE_VEC);
   MEM_STAT_REG(givs,TYPE_VEC);
   MEM_STAT_REG(givc,TYPE_VEC);
   
   R = m_resize(R,ip->k+1,ip->k);
   Q = m_resize(Q,ip->k,ip->b->dim);
   MEM_STAT_REG(R,TYPE_MAT);
   MEM_STAT_REG(Q,TYPE_MAT);		

   if (ip->x == VNULL) {  /* ip->x == 0 */
      ip->x = v_get(ip->b->dim);
      ip->shared_x = FALSE;
   }   

   v.dim = v.max_dim = ip->b->dim;      /* v and v1 are pointers to rows */
   v1.dim = v1.max_dim = ip->b->dim;  	/* of matrix Q */
   
   if (ip->Bx != (Fun_Ax)NULL) {    /* if precondition is defined */
      z = v_resize(z,ip->b->dim);
      MEM_STAT_REG(z,TYPE_VEC);
   }
   
   done = FALSE;
   for (ip->steps = 0; ip->steps < ip->limit; ) {

      /* restart */

      ip->Ax(ip->A_par,ip->x,u);    		/* u = A*x */
      v_sub(ip->b,u,u);		 		/* u = b - A*x */
      rr = u;				/* rr is a pointer only */
      
      if (ip->Bx) {
	 (ip->Bx)(ip->B_par,u,z);            /* tmp = B*(b-A*x)  */
	 rr = z;
      }
      
      nres = v_norm2(rr);
      if (ip->steps == 0) {
	 if (ip->info) ip->info(ip,nres,VNULL,VNULL);
	 ip->init_res = nres;
      }

      if ( nres == 0.0 ) {
	 done = TRUE;
	 break;
      }

      v.ve = Q->me[0];
      sv_mlt(1.0/nres,rr,&v);
      
      v_zero(r);
      v_zero(rhs);
      rhs->ve[0] = nres;

      for ( i = 0; i < ip->k && ip->steps < ip->limit; i++ ) {
	 ip->steps++;
	 v.ve = Q->me[i];	
	 (ip->Ax)(ip->A_par,&v,u);
	 rr = u;
	 if (ip->Bx) {
	    (ip->Bx)(ip->B_par,u,z);
	    rr = z;
	 }
	 
	 if (i < ip->k - 1) {
	    v1.ve = Q->me[i+1];
	    v_copy(rr,&v1);
	    for (j = 0; j <= i; j++) {
	       v.ve = Q->me[j];
	       /* r->ve[j] = in_prod(&v,rr); */
	       /* modified Gram-Schmidt algorithm */
	       r->ve[j] = in_prod(&v,&v1);
	       v_mltadd(&v1,&v,-r->ve[j],&v1);
	    }
	    
	    r->ve[i+1] = nres = v_norm2(&v1);
	    if (nres <= MACHEPS*ip->init_res) {
	       for (j = 0; j < i; j++) 
		 rot_vec(r,j,j+1,givc->ve[j],givs->ve[j],r);
	       set_col(R,i,r);
	       done = TRUE;
	       break;
	    }
	    sv_mlt(1.0/nres,&v1,&v1);
	 }
	 else {  /* i == ip->k - 1 */
	    /* Q->me[ip->k] need not be computed */

	    for (j = 0; j <= i; j++) {
	       v.ve = Q->me[j];
	       r->ve[j] = in_prod(&v,rr);
	    }
	    
	    nres = in_prod(rr,rr) - in_prod(r,r);
	    if (sqrt(fabs(nres)) <= MACHEPS*ip->init_res) { 
	       for (j = 0; j < i; j++) 
		 rot_vec(r,j,j+1,givc->ve[j],givs->ve[j],r);
	       set_col(R,i,r);
	       done = TRUE;
	       break;
	    }
	    if (nres < 0.0) { /* do restart */
	       i--; 
	       ip->steps--;
	       break;
	    } 
	    r->ve[i+1] = sqrt(nres);
	 }

	 /* QR update */

	 /* last_h = r->ve[i+1]; */ /* for test only */
	 for (j = 0; j < i; j++) 
	   rot_vec(r,j,j+1,givc->ve[j],givs->ve[j],r);
	 givens(r->ve[i],r->ve[i+1],&givc->ve[i],&givs->ve[i]);
	 rot_vec(r,i,i+1,givc->ve[i],givs->ve[i],r);
	 rot_vec(rhs,i,i+1,givc->ve[i],givs->ve[i],rhs);
	 
	 set_col(R,i,r);

	 nres = fabs((double) rhs->ve[i+1]);
	 if (ip->info) ip->info(ip,nres,VNULL,VNULL);
	 if ( ip->stop_crit(ip,nres,VNULL,VNULL) ) {
	    done = TRUE;
	    break;
	 }
      }
      
      /* use ixi submatrix of R */

      if (i >= ip->k) i = ip->k - 1;

      R = m_resize(R,i+1,i+1);
      rhs = v_resize(rhs,i+1);
      
      /* test only */
      /* test_gmres(ip,i,Q,R,givc,givs,last_h);  */
      
      Usolve(R,rhs,rhs,0.0); 	 /* solve a system: R*x = rhs */

      /* new approximation */

      for (j = 0; j <= i; j++) {
	 v.ve = Q->me[j]; 
	 v_mltadd(ip->x,&v,rhs->ve[j],ip->x);
      }

      if (done) break;

      /* back to old dimensions */

      rhs = v_resize(rhs,ip->k+1);
      R = m_resize(R,ip->k+1,ip->k);

   }

#ifdef THREADSAFE
   V_FREE(u);		V_FREE(r);	V_FREE(rhs);
   V_FREE(givs);	V_FREE(givc);	V_FREE(z);
   M_FREE(Q);		M_FREE(R);
#endif

   return ip->x;
}

/* iter_spgmres - a simple interface to iter_gmres */
#ifndef ANSI_C
VEC	*iter_spgmres(A,B,b,tol,x,k,limit,steps)
SPMAT	*A, *B;
VEC	*b, *x;
double	tol;
int *steps,k,limit;
#else
VEC	*iter_spgmres(SPMAT *A, SPMAT *B, VEC *b, double tol,
		      VEC *x, int k, 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->k = k;
   ip->limit = limit;
   ip->info = (Fun_info) NULL;
   ip->b = b;
   ip->eps = tol;
   ip->x = x;
   iter_gmres(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;		
}


/* for testing mgcr */
#ifndef ANSI_C
static void test_mgcr(ip,i,Q,R)
ITER *ip;
int i;
MAT *Q, *R;
#else
static void test_mgcr(ITER *ip, int i, MAT *Q, MAT *R)
#endif
{
   VEC vt, vt1;
   static MAT *R1=MNULL;
   static VEC *r=VNULL, *r1=VNULL;
   VEC *rr;
   int k,j;
   Real sm;
   
   
   /* check Q*Q^T = I */
   vt.dim = vt.max_dim = ip->b->dim;
   vt1.dim = vt1.max_dim = ip->b->dim;
   
   Q = m_resize(Q,i+1,ip->b->dim);
   R1 = m_resize(R1,i+1,i+1);
   r = v_resize(r,ip->b->dim);
   r1 = v_resize(r1,ip->b->dim);
   MEM_STAT_REG(R1,TYPE_MAT);
   MEM_STAT_REG(r,TYPE_VEC);
   MEM_STAT_REG(r1,TYPE_VEC);

   m_zero(R1);
   for (k=1; k <= i; k++)
     for (j=1; j <= i; j++) {
	vt.ve = Q->me[k];
	vt1.ve = Q->me[j];
	R1->me[k][j] = in_prod(&vt,&vt1);
     }
   for (j=1; j <= i; j++)