iternsym.c

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

     R1->me[j][j] -= 1.0;
#ifndef MEX
   if (m_norm_inf(R1) > MACHEPS*ip->b->dim)
     printf(" ! (mgcr:) m_norm_inf(Q*Q^T) = %g\n",m_norm_inf(R1));  
#endif

   /* check (r_i,Ap_j) = 0 for j <= i */
   
   ip->Ax(ip->A_par,ip->x,r);
   v_sub(ip->b,r,r);
   rr = r;
   if (ip->Bx) {
      ip->Bx(ip->B_par,r,r1);
      rr = r1;
   }
   
#ifndef MEX
   printf(" ||r|| = %g\n",v_norm2(rr));
#endif
   sm = 0.0;
   for (j = 1; j <= i; j++) {
      vt.ve = Q->me[j];
      sm = max(sm,in_prod(&vt,rr));
   }
#ifndef MEX
   if (sm >= MACHEPS*ip->b->dim)
     printf(" ! (mgcr:) max_j (r,Ap_j) = %g\n",sm);
#endif

}




/* 
  iter_mgcr -- modified generalized conjugate residual algorithm;
  fast version of GCR;
*/
#ifndef ANSI_C
VEC *iter_mgcr(ip)
ITER *ip;
#else
VEC *iter_mgcr(ITER *ip)
#endif
{
   STATIC VEC *As=VNULL, *beta=VNULL, *alpha=VNULL, *z=VNULL;
   STATIC MAT *N=MNULL, *H=MNULL;
   
   VEC *rr, v, s;  /* additional pointer and structures */
   Real nres;      /* norm of a residual */
   Real dd;        /* coefficient d_i */
   int i,j;
   int done;      /* if TRUE then stop the iterative process */
   int dim;       /* dimension of the problem */
   
   /* ip cannot be NULL */
   if (ip == INULL) error(E_NULL,"mgcr");
   /* Ax, b and stopping criterion must be given */
   if (! ip->Ax || ! ip->b || ! ip->stop_crit) 
     error(E_NULL,"mgcr");
   /* at least one direction vector must exist */
   if ( ip->k <= 0) error(E_BOUNDS,"mgcr");
   /* if the vector x is given then b and x must have the same dimension */
   if ( ip->x && ip->x->dim != ip->b->dim)
     error(E_SIZES,"mgcr");
   if (ip->eps <= 0.0) ip->eps = MACHEPS;
   
   dim = ip->b->dim;
   As = v_resize(As,dim);
   alpha = v_resize(alpha,ip->k);
   beta = v_resize(beta,ip->k);
   
   MEM_STAT_REG(As,TYPE_VEC);
   MEM_STAT_REG(alpha,TYPE_VEC);
   MEM_STAT_REG(beta,TYPE_VEC);
   
   H = m_resize(H,ip->k,ip->k);
   N = m_resize(N,ip->k,dim);
   
   MEM_STAT_REG(H,TYPE_MAT);
   MEM_STAT_REG(N,TYPE_MAT);
   
   /* if a preconditioner is defined */
   if (ip->Bx) {
      z = v_resize(z,dim);
      MEM_STAT_REG(z,TYPE_VEC);
   }
   
   /* if x is NULL then it is assumed that x has 
      entries with value zero */
   if ( ! ip->x ) {
      ip->x = v_get(ip->b->dim);
      ip->shared_x = FALSE;
   }
   
   /* v and s are additional pointers to rows of N */
   /* they must have the same dimension as rows of N */
   v.dim = v.max_dim = s.dim = s.max_dim = dim;
   
   
   done = FALSE;
   for (ip->steps = 0; ip->steps < ip->limit; ) {
      (*ip->Ax)(ip->A_par,ip->x,As);         /* As = A*x */
      v_sub(ip->b,As,As);                    /* As = b - A*x */
      rr = As;                               /* rr is an additional pointer */
      
      /* if a preconditioner is defined */
      if (ip->Bx) {
	 (*ip->Bx)(ip->B_par,As,z);               /* z = B*(b-A*x)  */
	 rr = z;                                  
      }
      
      /* norm of the residual */
      nres = v_norm2(rr);
      dd = nres;                            /* dd = ||r_i||  */
      
      /* check if the norm of the residual is zero */
      if (ip->steps == 0) {                
	 /* information for a user */
	 if (ip->info) (*ip->info)(ip,nres,As,rr); 
	 ip->init_res = fabs(nres);
      }

      if (nres == 0.0) { 
	 /* iterative process is finished */
	 done = TRUE; 
	 break;
      }
      
      /* save this residual in the first row of N */
      v.ve = N->me[0];
      v_copy(rr,&v);
      
      for (i = 0; i < ip->k && ip->steps < ip->limit; i++) {
	 ip->steps++;
	 v.ve = N->me[i];                /* pointer to a row of N (=s_i) */
	 /* note that we must use here &v, not v */
	 (*ip->Ax)(ip->A_par,&v,As); 
	 rr = As;                        /* As = A*s_i */
	 if (ip->Bx) {
	    (*ip->Bx)(ip->B_par,As,z);    /* z = B*A*s_i  */
	    rr = z;
	 }
	 
	 if (i < ip->k - 1) {
	    s.ve = N->me[i+1];         /* pointer to a row of N (=s_{i+1}) */
	    v_copy(rr,&s);                   /* s_{i+1} = B*A*s_i */
	    for (j = 0; j <= i-1; j++) {
	       v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
	       /* beta->ve[j] = in_prod(&v,rr); */      /* beta_{j,i} */
	       /* modified Gram-Schmidt algorithm */
	       beta->ve[j] = in_prod(&v,&s);  	         /* beta_{j,i} */
	                                 /* s_{i+1} -= beta_{j,i}*s_{j+1} */
	       v_mltadd(&s,&v,- beta->ve[j],&s);    
	    }
	    
	     /* beta_{i,i} = ||s_{i+1}||_2 */
	    beta->ve[i] = nres = v_norm2(&s);     
	    if ( nres <= MACHEPS*ip->init_res) { 
	       /* s_{i+1} == 0 */
	       i--;
	       done = TRUE;
	       break;
	    }
	    sv_mlt(1.0/nres,&s,&s);           /* normalize s_{i+1} */
	    
	    v.ve = N->me[0];
	    alpha->ve[i] = in_prod(&v,&s);     /* alpha_i = (s_0 , s_{i+1}) */
	    
	 }
	 else {
	    for (j = 0; j <= i-1; j++) {
	       v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
	       beta->ve[j] = in_prod(&v,rr);       /* beta_{j,i} */
	    }
	    
	    nres = in_prod(rr,rr);                 /* rr = B*A*s_{k-1} */
	    for (j = 0; j <= i-1; j++)
              nres -= beta->ve[j]*beta->ve[j];

	    if (sqrt(fabs(nres)) <= MACHEPS*ip->init_res)  {
	       /* s_k is zero */
	       i--;
	       done = TRUE;
	       break;
	    }
	    if (nres < 0.0) { /* do restart */
	       i--; 
	       ip->steps--;
	       break; 
	    }   
	    beta->ve[i] = sqrt(nres);         /* beta_{k-1,k-1} */
	    
	    v.ve = N->me[0];
	    alpha->ve[i] = in_prod(&v,rr); 
	    for (j = 0; j <= i-1; j++)
              alpha->ve[i] -= beta->ve[j]*alpha->ve[j];
	    alpha->ve[i] /= beta->ve[i];                /* alpha_{k-1} */
	    
	 }
	 
	 set_col(H,i,beta);

	 /* other method of computing dd */
	/* if (fabs((double)alpha->ve[i]) > dd)  {     
	    nres = - dd*dd + alpha->ve[i]*alpha->ve[i];
	    nres = sqrt((double) nres); 
	    if (ip->info) (*ip->info)(ip,-nres,VNULL,VNULL);  	
	    break;     
	 }  */
	 /* to avoid overflow/underflow in computing dd */
	 /* dd *= cos(asin((double)(alpha->ve[i]/dd))); */
	 
	 nres = alpha->ve[i]/dd;
	 if (fabs(nres-1.0) <= MACHEPS*ip->init_res) 
	   dd = 0.0;
	 else {
	    nres = 1.0 - nres*nres;
	    if (nres < 0.0) {
	       nres = sqrt((double) -nres); 
	       if (ip->info) (*ip->info)(ip,-dd*nres,VNULL,VNULL);  	
	       break;
	    }
	    dd *= sqrt((double) nres);  
	 }

	 if (ip->info) (*ip->info)(ip,dd,VNULL,VNULL);     
	 if ( ip->stop_crit(ip,dd,VNULL,VNULL) ) {
	    /* stopping criterion is satisfied */
	    done = TRUE;
	    break;
	 }
	 
      } /* end of for */
      
      if (i >= ip->k) i = ip->k - 1;
      
      /* use (i+1) by (i+1) submatrix of H */
      H = m_resize(H,i+1,i+1);
      alpha = v_resize(alpha,i+1);
      Usolve(H,alpha,alpha,0.0);       /* c_i is saved in alpha */
      
      for (j = 0; j <= i; j++) {
	 v.ve = N->me[j];
	 v_mltadd(ip->x,&v,alpha->ve[j],ip->x);
      }
      
      
      if (done) break;              /* stop the iterative process */
      alpha = v_resize(alpha,ip->k);
      H = m_resize(H,ip->k,ip->k);
      
   }  /* end of while */

#ifdef THREADSAFE
   V_FREE(As);		V_FREE(beta);	V_FREE(alpha);	V_FREE(z);
   M_FREE(N);		M_FREE(H);
#endif

   return ip->x;                    /* return the solution */
}



/* iter_spmgcr - a simple interface to iter_mgcr */
/* no preconditioner */
#ifndef ANSI_C
VEC	*iter_spmgcr(A,B,b,tol,x,k,limit,steps)
SPMAT	*A, *B;
VEC	*b, *x;
double	tol;
int *steps,k,limit;
#else
VEC	*iter_spmgcr(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_mgcr(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;		
}



/* 
  Conjugate gradients method for a normal equation
  a preconditioner B must be symmetric !!
*/
#ifndef ANSI_C
VEC  *iter_cgne(ip)
ITER *ip;
#else
VEC  *iter_cgne(ITER *ip)
#endif
{
   STATIC VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
   Real	alpha, beta, inner, old_inner, nres;
   VEC *rr1;   /* pointer only */
   
   if (ip == INULL)
     error(E_NULL,"iter_cgne");
   if (!ip->Ax || ! ip->ATx || !ip->b)
     error(E_NULL,"iter_cgne");
   if ( ip->x == ip->b )
     error(E_INSITU,"iter_cgne");
   if (!ip->stop_crit)
     error(E_NULL,"iter_cgne");
   
   if ( ip->eps <= 0.0 ) ip->eps = MACHEPS;
   
   r = v_resize(r,ip->b->dim);
   p = v_resize(p,ip->b->dim);
   q = v_resize(q,ip->b->dim);

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

   z = v_resize(z,ip->b->dim);
   MEM_STAT_REG(z,TYPE_VEC);

   if (ip->x) {
      if (ip->x->dim != ip->b->dim)
	error(E_SIZES,"iter_cgne");
      ip->Ax(ip->A_par,ip->x,p);    		/* p = A*x */
      v_sub(ip->b,p,z);		 		/* z = b - A*x */
   }
   else {  /* ip->x == 0 */
      ip->x = v_get(ip->b->dim);
      ip->shared_x = FALSE;
      v_copy(ip->b,z);
   }
   rr1 = z;
   if (ip->Bx) {
      (ip->Bx)(ip->B_par,rr1,p);
      rr1 = p;
   }
   (ip->ATx)(ip->AT_par,rr1,r);		/* r = A^T*B*(b-A*x)  */


   old_inner = 0.0;
   for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
   {
      rr1 = r;
      if ( ip->Bx ) {
	 (ip->Bx)(ip->B_par,r,z);		/* rr = B*r */
	 rr1 = z;
      }

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

      if ( ip->steps )	/* if ( ip->steps > 0 ) ... */
      {
	 beta = inner/old_inner;
	 p = v_mltadd(rr1,p,beta,p);
      }
      else		/* if ( ip->steps == 0 ) ... */
      {
	 beta = 0.0;
	 p = v_copy(rr1,p);
	 old_inner = 0.0;
      }
      (ip->Ax)(ip->A_par,p,q);     /* q = A*p */
      if (ip->Bx) {
	 (ip->Bx)(ip->B_par,q,z);
	 (ip->ATx)(ip->AT_par,z,q);
	 rr1 = q;			/* q = A^T*B*A*p */
      }
      else {
	 (ip->ATx)(ip->AT_par,q,z);	/* z = A^T*A*p */
	 rr1 = z;
      }

      alpha = inner/in_prod(rr1,p);
      v_mltadd(ip->x,p,alpha,ip->x);
      v_mltadd(r,rr1,-alpha,r);
      old_inner = inner;
   }

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

   return ip->x;
}

/* iter_spcgne -- a simple interface to iter_cgne() which 
   uses sparse matrix data structures
   -- assumes that B contains an actual preconditioner (or NULL)
   use always as follows:
      x = iter_spcgne(A,B,b,eps,x,limit,steps);
   or 
      x = iter_spcgne(A,B,b,eps,VNULL,limit,steps);
   In the second case the solution vector is created.
*/
#ifndef ANSI_C
VEC  *iter_spcgne(A,B,b,eps,x,limit,steps)
SPMAT	*A, *B;
VEC	*b, *x;
double	eps;
int *steps, limit;
#else
VEC  *iter_spcgne(SPMAT *A,SPMAT *B, VEC *b, double eps,
		  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;
   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->b = b;
   ip->eps = eps;
   ip->limit = limit;
   ip->x = x;
   iter_cgne(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;		
}