itertort.c

来自「C语言版本的矩阵库」· C语言 代码 · 共 692 行 · 第 1/2 页

   v_rand(v);
   v = iter_spcgs(An,NULL,yn,u,EPS,v,1000,&k);

   mem_stat_free(1);
   printf(" cgs: no. of iter.steps = %d\n",k);
   v_sub(v,xn,u);
   printf(" (cgs:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
	  v_norm2(u),EPS);
   


   /*** LSQR ***/

   notice("LSQR method (without preconditioning)");

   v_rand(u);
   v_free(ipns1->x);
   ipns1->x = u;
   ipns1->shared_x = TRUE;
   ipns1->info = NULL;
   mem_stat_mark(2);
   z = iter_lsqr(ipns1);
   
   v_sub(xn,z,v);
   k = ipns1->steps;
   printf(" lsqr: # of iter. steps = %d\n",k);
   printf(" (lsqr:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
	  v_norm2(v),EPS);

   v_rand(u);
   u = iter_splsqr(An,yn,EPS,u,1000,&k);
   mem_stat_free(2);
   
   v_sub(xn,u,v);
   printf(" splsqr: # of iter. steps = %d\n",k);
   printf(" (splsqr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n",	
	  v_norm2(v),EPS);



   /***** GMRES ********/

   notice("GMRES method with ICH preconditioning (nonsymmetric case)");

   v_zero(ipns->x);
/*   ipns->info = iter_std_info;  */
   ipns->info = NULL;  

   mem_stat_mark(2);
   z = iter_gmres(ipns);
   v_sub(xn,z,v);
   k = ipns->steps;
   printf(" gmres: # of iter. steps = %d\n",k);
   printf(" (gmres:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
	  v_norm2(v),EPS);

   notice("GMRES method without preconditioning (nonsymmetric case)");
   V_FREE(v);
   v = iter_spgmres(An,NULL,yn,EPS,VNULL,10,1004,&k);
   mem_stat_free(2);
   
   v_sub(xn,v,v);
   printf(" spgmres: # of iter. steps = %d\n",k);
   printf(" (spgmres:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
	  v_norm2(v),EPS);



   /**** MGCR *****/

   notice("MGCR method with ICH preconditioning (nonsymmetric case)");

   v_zero(ipns->x);
   mem_stat_mark(2);
   z = iter_mgcr(ipns);
   v_sub(xn,z,v);
   k = ipns->steps;
   printf(" mgcr: # of iter. steps = %d\n",k);
   printf(" (mgcr:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
	  v_norm2(v),EPS);

   notice("MGCR method without  preconditioning (nonsymmetric case)");
   V_FREE(v);
   v = iter_spmgcr(An,NULL,yn,EPS,VNULL,10,1004,&k);
   mem_stat_free(2);
   
   v_sub(xn,v,v);
   printf(" spmgcr: # of iter. steps = %d\n",k);
   printf(" (spmgcr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n",
	  v_norm2(v),EPS);


   /***** ARNOLDI METHOD ********/


   notice("arnoldi method");

   kk = ipns1->k = KK;
   Q = m_get(kk,x->dim);
   Q1 = m_get(kk,x->dim);
   H = m_get(kk,kk);
   v_rand(u);
   ipns1->x = u;
   ipns1->shared_x = TRUE;
   mem_stat_mark(3);
   iter_arnoldi_iref(ipns1,&hh,Q,H);
   mem_stat_free(3);

   /* check the equality:
      Q*A*Q^T = H; */

   vt.dim = vt.max_dim = x->dim;
   vt1.dim = vt1.max_dim = x->dim;
   for (j=0; j < kk; j++) {
      vt.ve = Q->me[j];
      vt1.ve = Q1->me[j];
      sp_mv_mlt(An,&vt,&vt1);
   }
   H1 = m_get(kk,kk);
   mmtr_mlt(Q,Q1,H1);
   m_sub(H,H1,H1);
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (arnoldi_iref) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);

   /* check Q*Q^T = I  */

   mmtr_mlt(Q,Q,H1);
   for (j=0; j < kk; j++)
     H1->me[j][j] -= 1.0;
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (arnoldi_iref) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);

   ipns1->x = u;
   ipns1->shared_x = TRUE;
   mem_stat_mark(3);
   iter_arnoldi(ipns1,&hh,Q,H);
   mem_stat_free(3);

   /* check the equality:
      Q*A*Q^T = H; */

   vt.dim = vt.max_dim = x->dim;
   vt1.dim = vt1.max_dim = x->dim;
   for (j=0; j < kk; j++) {
      vt.ve = Q->me[j];
      vt1.ve = Q1->me[j];
      sp_mv_mlt(An,&vt,&vt1);
   }

   mmtr_mlt(Q,Q1,H1);
   m_sub(H,H1,H1);
  if (m_norm_inf(H1) > MACHEPS*x->dim)  
     printf(" (arnoldi) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);
   /* check Q*Q^T = I  */
   mmtr_mlt(Q,Q,H1);
   for (j=0; j < kk; j++)
     H1->me[j][j] -= 1.0;
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (arnoldi) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);

   v_rand(u);
   mem_stat_mark(3);
   iter_sparnoldi(An,u,kk,&hh,Q,H);
   mem_stat_free(3);

   /* check the equality:
      Q*A*Q^T = H; */

   vt.dim = vt.max_dim = x->dim;
   vt1.dim = vt1.max_dim = x->dim;
   for (j=0; j < kk; j++) {
      vt.ve = Q->me[j];
      vt1.ve = Q1->me[j];
      sp_mv_mlt(An,&vt,&vt1);
   }

   mmtr_mlt(Q,Q1,H1);
   m_sub(H,H1,H1);
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (sparnoldi) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);
   /* check Q*Q^T = I  */
   mmtr_mlt(Q,Q,H1);
   for (j=0; j < kk; j++)
     H1->me[j][j] -= 1.0;
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (sparnoldi) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);



   /****** LANCZOS METHOD ******/

   notice("lanczos method");
   kk = ipns1->k; 
   Q = m_resize(Q,kk,x->dim);
   Q1 = m_resize(Q1,kk,x->dim);
   H = m_resize(H,kk,kk);
   ips1->k = kk;
   v_rand(u);
   v_free(ips1->x);
   ips1->x = u;
   ips1->shared_x = TRUE;

   mem_stat_mark(3);
   iter_lanczos(ips1,x,y,&hh,Q);
   mem_stat_free(3);

   /* check the equality:
      Q*A*Q^T = H; */

   vt.dim = vt1.dim = Q->n;
   vt.max_dim = vt1.max_dim = Q->max_n;
   Q1 = m_resize(Q1,Q->m,Q->n);
   for (j=0; j < Q->m; j++) {
      vt.ve = Q->me[j];
      vt1.ve = Q1->me[j];
      sp_mv_mlt(A,&vt,&vt1);
   }
   H1 = m_resize(H1,Q->m,Q->m);
   H = m_resize(H,Q->m,Q->m);
   mmtr_mlt(Q,Q1,H1);

   m_zero(H);
   for (j=0; j < Q->m-1; j++) {
      H->me[j][j] = x->ve[j];
      H->me[j][j+1] = H->me[j+1][j] = y->ve[j];
   }
   H->me[Q->m-1][Q->m-1] = x->ve[Q->m-1];

   m_sub(H,H1,H1);
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (lanczos) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);

   /* check Q*Q^T = I  */

   mmtr_mlt(Q,Q,H1);
   for (j=0; j < Q->m; j++)
     H1->me[j][j] -= 1.0;
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (lanczos) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);

   mem_stat_mark(3);
   v_rand(u);
   iter_splanczos(A,kk,u,x,y,&hh,Q);
   mem_stat_free(3);

   /* check the equality:
      Q*A*Q^T = H; */

   vt.dim = vt1.dim = Q->n;
   vt.max_dim = vt1.max_dim = Q->max_n;
   Q1 = m_resize(Q1,Q->m,Q->n);
   for (j=0; j < Q->m; j++) {
      vt.ve = Q->me[j];
      vt1.ve = Q1->me[j];
      sp_mv_mlt(A,&vt,&vt1);
   }
   H1 = m_resize(H1,Q->m,Q->m);
   H = m_resize(H,Q->m,Q->m);
   mmtr_mlt(Q,Q1,H1);
   for (j=0; j < Q->m-1; j++) {
      H->me[j][j] = x->ve[j];
      H->me[j][j+1] = H->me[j+1][j] = y->ve[j];
   }
   H->me[Q->m-1][Q->m-1] = x->ve[Q->m-1];

   m_sub(H,H1,H1);
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (splanczos) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);
   /* check Q*Q^T = I  */
   mmtr_mlt(Q,Q,H1);
   for (j=0; j < Q->m; j++)
     H1->me[j][j] -= 1.0;
   if (m_norm_inf(H1) > MACHEPS*x->dim)
     printf(" (splanczos) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
	    m_norm_inf(H1),MACHEPS);



   /***** LANCZOS2 ****/

   notice("lanczos2 method");
   kk = 50;  		/* # of dir. vectors */
   ips1->k = kk;
   v_rand(u);
   ips1->x = u;
   ips1->shared_x = TRUE;

   for ( i = 0; i < xn->dim; i++ )
	xn->ve[i] = i;
   iter_Ax(ips1,Dv_mlt,xn);
   mem_stat_mark(3);
   iter_lanczos2(ips1,y,v);
   mem_stat_free(3);

   printf("# Number of steps of Lanczos algorithm = %d\n", kk);
   printf("# Exact eigenvalues are 0, 1, 2, ..., %d\n",ANON-1);
   printf("# Extreme eigenvalues should be accurate; \n");
   printf("# interior values usually are not.\n");
   printf("# approx e-vals =\n");	v_output(y);
   printf("# Error in estimate of bottom e-vec (Lanczos) = %g\n",
	  fabs(v->ve[0]));

   mem_stat_mark(3);
   v_rand(u);
   iter_splanczos2(A,kk,u,y,v);
   mem_stat_free(3);


   /***** FINISHING *******/

   notice("release ITER variables");
   
   M_FREE(Q);
   M_FREE(Q1);
   M_FREE(H);
   M_FREE(H1);

   ITER_FREE(ipns);
   ITER_FREE(ips);
   ITER_FREE(ipns1);
   ITER_FREE(ips1);
   SP_FREE(A);
   SP_FREE(B);
   SP_FREE(An);
   SP_FREE(Bn);
   
   V_FREE(x);
   V_FREE(y);
   V_FREE(u);
   V_FREE(v); 
   V_FREE(xn);
   V_FREE(yn);

   printf("# Done testing (%s)\n",argv[0]);
   mem_info();
}