home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / itertort < prev    next >
Encoding:
Text File  |  1994-01-13  |  15.8 KB  |  692 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*  iter_tort.c  16/09/93 */
  28.  
  29. /*
  30.   This file contains tests for the iterative part of Meschach
  31. */
  32.  
  33. #include    <stdio.h>
  34. #include    <math.h>
  35. #include    "matrix2.h"
  36. #include    "sparse2.h"
  37. #include    "iter.h"
  38.  
  39. #define    errmesg(mesg)    printf("Error: %s error: line %d\n",mesg,__LINE__)
  40. #define notice(mesg)    printf("# Testing %s...\n",mesg);
  41.   
  42.   /* for iterative methods */
  43.   
  44. #if REAL == DOUBLE
  45. #define    EPS    1e-7
  46. #define KK    20
  47. #elif REAL == FLOAT
  48. #define EPS   1e-5
  49. #define KK    8
  50. #endif
  51.  
  52. #define ANON  513
  53. #define ASYM  ANON   
  54.  
  55.   
  56. static VEC *ex_sol = VNULL;
  57.  
  58. /* new iter information */
  59. void iter_mod_info(ip,nres,res,Bres)
  60. ITER *ip;
  61. double nres;
  62. VEC *res, *Bres;
  63. {
  64.    static VEC *tmp;
  65.  
  66.    if (ip->b == VNULL) return;
  67.    tmp = v_resize(tmp,ip->b->dim);
  68.    MEM_STAT_REG(tmp,TYPE_VEC);
  69.  
  70.    if (nres >= 0.0) {
  71.       printf(" %d. residual = %g\n",ip->steps,nres);
  72.    }
  73.    else 
  74.      printf(" %d. residual = %g (WARNING !!! should be >= 0) \n",
  75.         ip->steps,nres);
  76.    if (ex_sol != VNULL)
  77.      printf("    ||u_ex - u_approx||_2 = %g\n",
  78.         v_norm2(v_sub(ip->x,ex_sol,tmp)));
  79. }
  80.  
  81.  
  82. /* out = A^T*A*x */
  83. VEC *norm_equ(A,x,out)
  84. SPMAT *A;
  85. VEC *x, *out;
  86. {
  87.    static VEC * tmp;
  88.  
  89.    tmp = v_resize(tmp,x->dim);
  90.    MEM_STAT_REG(tmp,TYPE_VEC);
  91.    sp_mv_mlt(A,x,tmp);
  92.    sp_vm_mlt(A,tmp,out);
  93.    return out;
  94.  
  95. }
  96.  
  97.  
  98. /* 
  99.   make symmetric preconditioner for nonsymmetric matrix A;
  100.    B = 0.5*(A+A^T) and then B is factorized using 
  101.    incomplete Choleski factorization
  102. */
  103.  
  104. SPMAT *gen_sym_precond(A)
  105. SPMAT *A;
  106. {
  107.    SPMAT *B;
  108.    SPROW *row;
  109.    int i,j,k;
  110.    Real val;
  111.    
  112.    B = sp_get(A->m,A->n,A->row[0].maxlen);
  113.    for (i=0; i < A->m; i++) {
  114.       row = &(A->row[i]);
  115.       for (j = 0; j < row->len; j++) {
  116.     k = row->elt[j].col;
  117.     if (i != k) {
  118.        val = 0.5*(sp_get_val(A,i,k) + sp_get_val(A,k,i));
  119.        sp_set_val(B,i,k,val);
  120.        sp_set_val(B,k,i,val);
  121.     }
  122.     else { /* i == k */
  123.       val = sp_get_val(A,i,i);
  124.       sp_set_val(B,i,i,val);
  125.        }
  126.      }
  127.    }
  128.  
  129.    spICHfactor(B);
  130.    return B;
  131. }
  132.  
  133. /* Dv_mlt -- diagonal by vector multiply; the diagonal matrix is represented
  134.         by a vector d */
  135. VEC    *Dv_mlt(d, x, out)
  136. VEC    *d, *x, *out;
  137. {
  138.     int        i;
  139.  
  140.     if ( ! d || ! x )
  141.     error(E_NULL,"Dv_mlt");
  142.     if ( d->dim != x->dim )
  143.     error(E_SIZES,"Dv_mlt");
  144.     out = v_resize(out,x->dim);
  145.  
  146.     for ( i = 0; i < x->dim; i++ )
  147.     out->ve[i] = d->ve[i]*x->ve[i];
  148.  
  149.     return out;
  150. }
  151.  
  152.  
  153.  
  154. /************************************************/
  155. void    main(argc, argv)
  156. int    argc;
  157. char    *argv[];
  158. {
  159.    VEC        *x, *y, *z, *u, *v, *xn, *yn;
  160.    SPMAT    *A = NULL, *B = NULL;
  161.    SPMAT    *An = NULL, *Bn = NULL;
  162.    int        i, k, kk, j;
  163.    ITER        *ips, *ips1, *ipns, *ipns1;
  164.    MAT         *Q, *H, *Q1, *H1;
  165.    VEC         vt, vt1;
  166.    Real        hh;
  167.  
  168.  
  169.    mem_info_on(TRUE);
  170.    notice("allocating sparse matrices");
  171.    
  172.    printf(" dim of A = %dx%d\n",ASYM,ASYM);
  173.    
  174.    A = iter_gen_sym(ASYM,8);   
  175.    B = sp_copy(A);
  176.    spICHfactor(B);
  177.    
  178.    u = v_get(A->n);
  179.    x = v_get(A->n);
  180.    y = v_get(A->n);
  181.    v = v_get(A->n);
  182.  
  183.    v_rand(x);
  184.    sp_mv_mlt(A,x,y);
  185.    ex_sol = x;
  186.    
  187.    notice(" initialize ITER variables");
  188.    /* ips for symmetric matrices with precondition */
  189.    ips = iter_get(A->m,A->n);
  190.  
  191.    /*  printf(" ips:\n");
  192.    iter_dump(stdout,ips);   */
  193.  
  194.    ips->limit = 500;
  195.    ips->eps = EPS;
  196.    
  197.    iter_Ax(ips,sp_mv_mlt,A);
  198.    iter_Bx(ips,spCHsolve,B);
  199.  
  200.    ips->b = v_copy(y,ips->b);
  201.    v_rand(ips->x);
  202.    /* test of iter_resize */
  203.    ips = iter_resize(ips,2*A->m,2*A->n);
  204.    ips = iter_resize(ips,A->m,A->n);
  205.  
  206.    /*  printf(" ips:\n");
  207.    iter_dump(stdout,ips); */
  208.    
  209.    /* ips1 for symmetric matrices without precondition */
  210.    ips1 = iter_get(0,0);
  211.    /*   printf(" ips1:\n");
  212.    iter_dump(stdout,ips1);   */
  213.    ITER_FREE(ips1);
  214.  
  215.    ips1 = iter_copy2(ips,ips1);
  216.    iter_Bx(ips1,NULL,NULL);
  217.    ips1->b = ips->b;
  218.    ips1->shared_b = TRUE;
  219.    /*    printf(" ips1:\n");
  220.    iter_dump(stdout,ips1);   */
  221.  
  222.    /* ipns for nonsymetric matrices with precondition */
  223.    ipns = iter_copy(ips,INULL);
  224.    ipns->k = KK;
  225.    ipns->limit = 500;
  226.    ipns->info = NULL;
  227.  
  228.    An = iter_gen_nonsym_posdef(ANON,8);   
  229.    Bn = gen_sym_precond(An);
  230.    xn = v_get(An->n);
  231.    yn = v_get(An->n);
  232.    v_rand(xn);
  233.    sp_mv_mlt(An,xn,yn);
  234.    ipns->b = v_copy(yn,ipns->b);
  235.  
  236.    iter_Ax(ipns, sp_mv_mlt,An);
  237.    iter_ATx(ipns, sp_vm_mlt,An);
  238.    iter_Bx(ipns, spCHsolve,Bn);
  239.  
  240.    /*  printf(" ipns:\n");
  241.    iter_dump(stdout,ipns); */
  242.    
  243.    /* ipns1 for nonsymmetric matrices without precondition */
  244.    ipns1 = iter_copy2(ipns,INULL);
  245.    ipns1->b = ipns->b;
  246.    ipns1->shared_b = TRUE;
  247.    iter_Bx(ipns1,NULL,NULL);
  248.  
  249.    /*   printf(" ipns1:\n");
  250.    iter_dump(stdout,ipns1);  */
  251.  
  252.  
  253.    /*******  CG  ********/
  254.  
  255.    notice(" CG method without preconditioning");
  256.    ips1->info = NULL;
  257.    mem_stat_mark(1);
  258.    iter_cg(ips1);
  259.  
  260.    k = ips1->steps;
  261.    z = ips1->x;
  262.    printf(" cg: no. of iter.steps = %d\n",k);
  263.    v_sub(z,x,u);
  264.    printf(" (cg:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n",
  265.       v_norm2(u),EPS);
  266.    
  267.    notice(" CG method with ICH preconditioning");
  268.  
  269.    ips->info = NULL;
  270.    v_zero(ips->x);  
  271.    iter_cg(ips);  
  272.  
  273.    k = ips->steps;
  274.    printf(" cg: no. of iter.steps = %d\n",k);
  275.    v_sub(ips->x,x,u);
  276.    printf(" (cg:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n",
  277.       v_norm2(u),EPS);
  278.    
  279.    V_FREE(v);
  280.    if ((v = iter_spcg(A,B,y,EPS,VNULL,1000,&k)) == VNULL)
  281.      errmesg("CG method with precond.: NULL solution"); 
  282.    
  283.    v_sub(ips->x,v,u);
  284.    if (v_norm2(u) >= EPS) {
  285.       errmesg("CG method with precond.: different solutions");
  286.       printf(" diff. = %g\n",v_norm2(u));
  287.    }   
  288.    
  289.  
  290.    mem_stat_free(1);
  291.    printf(" spcg: # of iter. steps = %d\n",k);
  292.    v_sub(v,x,u);
  293.    printf(" (spcg:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  294.       v_norm2(u),EPS);  
  295.  
  296.  
  297.    /***  CG FOR NORMAL EQUATION *****/
  298.  
  299.    notice("CGNE method with ICH preconditioning (nonsymmetric case)");
  300.  
  301.    /* ipns->info = iter_std_info;  */
  302.    ipns->info = NULL;
  303.    v_zero(ipns->x);
  304.  
  305.    mem_stat_mark(1);
  306.    iter_cgne(ipns);
  307.  
  308.    k = ipns->steps;
  309.    z = ipns->x;
  310.    printf(" cgne: no. of iter.steps = %d\n",k);
  311.    v_sub(z,xn,u);
  312.    printf(" (cgne:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  313.       v_norm2(u),EPS);
  314.  
  315.    notice("CGNE method without preconditioning (nonsymmetric case)");
  316.  
  317.    v_rand(u);
  318.    u = iter_spcgne(An,NULL,yn,EPS,u,1000,&k);
  319.  
  320.    mem_stat_free(1);
  321.    printf(" spcgne: no. of iter.steps = %d\n",k);
  322.    v_sub(u,xn,u);
  323.    printf(" (spcgne:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  324.       v_norm2(u),EPS);
  325.  
  326.    /***  CGS  *****/
  327.  
  328.    notice("CGS method with ICH preconditioning (nonsymmetric case)");
  329.  
  330.    v_zero(ipns->x);   /* new init guess == 0 */
  331.  
  332.    mem_stat_mark(1);
  333.    ipns->info = NULL;
  334.    v_rand(u);
  335.    iter_cgs(ipns,u);
  336.  
  337.    k = ipns->steps;
  338.    z = ipns->x;
  339.    printf(" cgs: no. of iter.steps = %d\n",k);
  340.    v_sub(z,xn,u);
  341.    printf(" (cgs:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  342.       v_norm2(u),EPS);
  343.  
  344.    notice("CGS method without preconditioning (nonsymmetric case)");
  345.  
  346.    v_rand(u);
  347.    v_rand(v);
  348.    v = iter_spcgs(An,NULL,yn,u,EPS,v,1000,&k);
  349.  
  350.    mem_stat_free(1);
  351.    printf(" cgs: no. of iter.steps = %d\n",k);
  352.    v_sub(v,xn,u);
  353.    printf(" (cgs:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  354.       v_norm2(u),EPS);
  355.    
  356.  
  357.  
  358.    /*** LSQR ***/
  359.  
  360.    notice("LSQR method (without preconditioning)");
  361.  
  362.    v_rand(u);
  363.    v_free(ipns1->x);
  364.    ipns1->x = u;
  365.    ipns1->shared_x = TRUE;
  366.    ipns1->info = NULL;
  367.    mem_stat_mark(2);
  368.    z = iter_lsqr(ipns1);
  369.    
  370.    v_sub(xn,z,v);
  371.    k = ipns1->steps;
  372.    printf(" lsqr: # of iter. steps = %d\n",k);
  373.    printf(" (lsqr:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  374.       v_norm2(v),EPS);
  375.  
  376.    v_rand(u);
  377.    u = iter_splsqr(An,yn,EPS,u,1000,&k);
  378.    mem_stat_free(2);
  379.    
  380.    v_sub(xn,u,v);
  381.    printf(" splsqr: # of iter. steps = %d\n",k);
  382.    printf(" (splsqr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n",    
  383.       v_norm2(v),EPS);
  384.  
  385.  
  386.  
  387.    /***** GMRES ********/
  388.  
  389.    notice("GMRES method with ICH preconditioning (nonsymmetric case)");
  390.  
  391.    v_zero(ipns->x);
  392. /*   ipns->info = iter_std_info;  */
  393.    ipns->info = NULL;  
  394.  
  395.    mem_stat_mark(2);
  396.    z = iter_gmres(ipns);
  397.    v_sub(xn,z,v);
  398.    k = ipns->steps;
  399.    printf(" gmres: # of iter. steps = %d\n",k);
  400.    printf(" (gmres:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  401.       v_norm2(v),EPS);
  402.  
  403.    notice("GMRES method without preconditioning (nonsymmetric case)");
  404.    V_FREE(v);
  405.    v = iter_spgmres(An,NULL,yn,EPS,VNULL,10,1004,&k);
  406.    mem_stat_free(2);
  407.    
  408.    v_sub(xn,v,v);
  409.    printf(" spgmres: # of iter. steps = %d\n",k);
  410.    printf(" (spgmres:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  411.       v_norm2(v),EPS);
  412.  
  413.  
  414.  
  415.    /**** MGCR *****/
  416.  
  417.    notice("MGCR method with ICH preconditioning (nonsymmetric case)");
  418.  
  419.    v_zero(ipns->x);
  420.    mem_stat_mark(2);
  421.    z = iter_mgcr(ipns);
  422.    v_sub(xn,z,v);
  423.    k = ipns->steps;
  424.    printf(" mgcr: # of iter. steps = %d\n",k);
  425.    printf(" (mgcr:) ||u_ex - u_approx||_2 = %g  [EPS = %g]\n",
  426.       v_norm2(v),EPS);
  427.  
  428.    notice("MGCR method without  preconditioning (nonsymmetric case)");
  429.    V_FREE(v);
  430.    v = iter_spmgcr(An,NULL,yn,EPS,VNULL,10,1004,&k);
  431.    mem_stat_free(2);
  432.    
  433.    v_sub(xn,v,v);
  434.    printf(" spmgcr: # of iter. steps = %d\n",k);
  435.    printf(" (spmgcr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n",
  436.       v_norm2(v),EPS);
  437.  
  438.  
  439.    /***** ARNOLDI METHOD ********/
  440.  
  441.  
  442.    notice("arnoldi method");
  443.  
  444.    kk = ipns1->k = KK;
  445.    Q = m_get(kk,x->dim);
  446.    Q1 = m_get(kk,x->dim);
  447.    H = m_get(kk,kk);
  448.    v_rand(u);
  449.    ipns1->x = u;
  450.    ipns1->shared_x = TRUE;
  451.    mem_stat_mark(3);
  452.    iter_arnoldi_iref(ipns1,&hh,Q,H);
  453.    mem_stat_free(3);
  454.  
  455.    /* check the equality:
  456.       Q*A*Q^T = H; */
  457.  
  458.    vt.dim = vt.max_dim = x->dim;
  459.    vt1.dim = vt1.max_dim = x->dim;
  460.    for (j=0; j < kk; j++) {
  461.       vt.ve = Q->me[j];
  462.       vt1.ve = Q1->me[j];
  463.       sp_mv_mlt(An,&vt,&vt1);
  464.    }
  465.    H1 = m_get(kk,kk);
  466.    mmtr_mlt(Q,Q1,H1);
  467.    m_sub(H,H1,H1);
  468.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  469.      printf(" (arnoldi_iref) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
  470.         m_norm_inf(H1),MACHEPS);
  471.  
  472.    /* check Q*Q^T = I  */
  473.  
  474.    mmtr_mlt(Q,Q,H1);
  475.    for (j=0; j < kk; j++)
  476.      H1->me[j][j] -= 1.0;
  477.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  478.      printf(" (arnoldi_iref) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
  479.         m_norm_inf(H1),MACHEPS);
  480.  
  481.    ipns1->x = u;
  482.    ipns1->shared_x = TRUE;
  483.    mem_stat_mark(3);
  484.    iter_arnoldi(ipns1,&hh,Q,H);
  485.    mem_stat_free(3);
  486.  
  487.    /* check the equality:
  488.       Q*A*Q^T = H; */
  489.  
  490.    vt.dim = vt.max_dim = x->dim;
  491.    vt1.dim = vt1.max_dim = x->dim;
  492.    for (j=0; j < kk; j++) {
  493.       vt.ve = Q->me[j];
  494.       vt1.ve = Q1->me[j];
  495.       sp_mv_mlt(An,&vt,&vt1);
  496.    }
  497.  
  498.    mmtr_mlt(Q,Q1,H1);
  499.    m_sub(H,H1,H1);
  500.   if (m_norm_inf(H1) > MACHEPS*x->dim)  
  501.      printf(" (arnoldi) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
  502.         m_norm_inf(H1),MACHEPS);
  503.    /* check Q*Q^T = I  */
  504.    mmtr_mlt(Q,Q,H1);
  505.    for (j=0; j < kk; j++)
  506.      H1->me[j][j] -= 1.0;
  507.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  508.      printf(" (arnoldi) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
  509.         m_norm_inf(H1),MACHEPS);
  510.  
  511.    v_rand(u);
  512.    mem_stat_mark(3);
  513.    iter_sparnoldi(An,u,kk,&hh,Q,H);
  514.    mem_stat_free(3);
  515.  
  516.    /* check the equality:
  517.       Q*A*Q^T = H; */
  518.  
  519.    vt.dim = vt.max_dim = x->dim;
  520.    vt1.dim = vt1.max_dim = x->dim;
  521.    for (j=0; j < kk; j++) {
  522.       vt.ve = Q->me[j];
  523.       vt1.ve = Q1->me[j];
  524.       sp_mv_mlt(An,&vt,&vt1);
  525.    }
  526.  
  527.    mmtr_mlt(Q,Q1,H1);
  528.    m_sub(H,H1,H1);
  529.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  530.      printf(" (sparnoldi) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
  531.         m_norm_inf(H1),MACHEPS);
  532.    /* check Q*Q^T = I  */
  533.    mmtr_mlt(Q,Q,H1);
  534.    for (j=0; j < kk; j++)
  535.      H1->me[j][j] -= 1.0;
  536.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  537.      printf(" (sparnoldi) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
  538.         m_norm_inf(H1),MACHEPS);
  539.  
  540.  
  541.  
  542.    /****** LANCZOS METHOD ******/
  543.  
  544.    notice("lanczos method");
  545.    kk = ipns1->k; 
  546.    Q = m_resize(Q,kk,x->dim);
  547.    Q1 = m_resize(Q1,kk,x->dim);
  548.    H = m_resize(H,kk,kk);
  549.    ips1->k = kk;
  550.    v_rand(u);
  551.    v_free(ips1->x);
  552.    ips1->x = u;
  553.    ips1->shared_x = TRUE;
  554.  
  555.    mem_stat_mark(3);
  556.    iter_lanczos(ips1,x,y,&hh,Q);
  557.    mem_stat_free(3);
  558.  
  559.    /* check the equality:
  560.       Q*A*Q^T = H; */
  561.  
  562.    vt.dim = vt1.dim = Q->n;
  563.    vt.max_dim = vt1.max_dim = Q->max_n;
  564.    Q1 = m_resize(Q1,Q->m,Q->n);
  565.    for (j=0; j < Q->m; j++) {
  566.       vt.ve = Q->me[j];
  567.       vt1.ve = Q1->me[j];
  568.       sp_mv_mlt(A,&vt,&vt1);
  569.    }
  570.    H1 = m_resize(H1,Q->m,Q->m);
  571.    H = m_resize(H,Q->m,Q->m);
  572.    mmtr_mlt(Q,Q1,H1);
  573.  
  574.    m_zero(H);
  575.    for (j=0; j < Q->m-1; j++) {
  576.       H->me[j][j] = x->ve[j];
  577.       H->me[j][j+1] = H->me[j+1][j] = y->ve[j];
  578.    }
  579.    H->me[Q->m-1][Q->m-1] = x->ve[Q->m-1];
  580.  
  581.    m_sub(H,H1,H1);
  582.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  583.      printf(" (lanczos) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
  584.         m_norm_inf(H1),MACHEPS);
  585.  
  586.    /* check Q*Q^T = I  */
  587.  
  588.    mmtr_mlt(Q,Q,H1);
  589.    for (j=0; j < Q->m; j++)
  590.      H1->me[j][j] -= 1.0;
  591.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  592.      printf(" (lanczos) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
  593.         m_norm_inf(H1),MACHEPS);
  594.  
  595.    mem_stat_mark(3);
  596.    v_rand(u);
  597.    iter_splanczos(A,kk,u,x,y,&hh,Q);
  598.    mem_stat_free(3);
  599.  
  600.    /* check the equality:
  601.       Q*A*Q^T = H; */
  602.  
  603.    vt.dim = vt1.dim = Q->n;
  604.    vt.max_dim = vt1.max_dim = Q->max_n;
  605.    Q1 = m_resize(Q1,Q->m,Q->n);
  606.    for (j=0; j < Q->m; j++) {
  607.       vt.ve = Q->me[j];
  608.       vt1.ve = Q1->me[j];
  609.       sp_mv_mlt(A,&vt,&vt1);
  610.    }
  611.    H1 = m_resize(H1,Q->m,Q->m);
  612.    H = m_resize(H,Q->m,Q->m);
  613.    mmtr_mlt(Q,Q1,H1);
  614.    for (j=0; j < Q->m-1; j++) {
  615.       H->me[j][j] = x->ve[j];
  616.       H->me[j][j+1] = H->me[j+1][j] = y->ve[j];
  617.    }
  618.    H->me[Q->m-1][Q->m-1] = x->ve[Q->m-1];
  619.  
  620.    m_sub(H,H1,H1);
  621.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  622.      printf(" (splanczos) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
  623.         m_norm_inf(H1),MACHEPS);
  624.    /* check Q*Q^T = I  */
  625.    mmtr_mlt(Q,Q,H1);
  626.    for (j=0; j < Q->m; j++)
  627.      H1->me[j][j] -= 1.0;
  628.    if (m_norm_inf(H1) > MACHEPS*x->dim)
  629.      printf(" (splanczos) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
  630.         m_norm_inf(H1),MACHEPS);
  631.  
  632.  
  633.  
  634.    /***** LANCZOS2 ****/
  635.  
  636.    notice("lanczos2 method");
  637.    kk = 50;          /* # of dir. vectors */
  638.    ips1->k = kk;
  639.    v_rand(u);
  640.    ips1->x = u;
  641.    ips1->shared_x = TRUE;
  642.  
  643.    for ( i = 0; i < xn->dim; i++ )
  644.     xn->ve[i] = i;
  645.    iter_Ax(ips1,Dv_mlt,xn);
  646.    mem_stat_mark(3);
  647.    iter_lanczos2(ips1,y,v);
  648.    mem_stat_free(3);
  649.  
  650.    printf("# Number of steps of Lanczos algorithm = %d\n", kk);
  651.    printf("# Exact eigenvalues are 0, 1, 2, ..., %d\n",ANON-1);
  652.    printf("# Extreme eigenvalues should be accurate; \n");
  653.    printf("# interior values usually are not.\n");
  654.    printf("# approx e-vals =\n");    v_output(y);
  655.    printf("# Error in estimate of bottom e-vec (Lanczos) = %g\n",
  656.       fabs(v->ve[0]));
  657.  
  658.    mem_stat_mark(3);
  659.    v_rand(u);
  660.    iter_splanczos2(A,kk,u,y,v);
  661.    mem_stat_free(3);
  662.  
  663.  
  664.    /***** FINISHING *******/
  665.  
  666.    notice("release ITER variables");
  667.    
  668.    M_FREE(Q);
  669.    M_FREE(Q1);
  670.    M_FREE(H);
  671.    M_FREE(H1);
  672.  
  673.    ITER_FREE(ipns);
  674.    ITER_FREE(ips);
  675.    ITER_FREE(ipns1);
  676.    ITER_FREE(ips1);
  677.    SP_FREE(A);
  678.    SP_FREE(B);
  679.    SP_FREE(An);
  680.    SP_FREE(Bn);
  681.    
  682.    V_FREE(x);
  683.    V_FREE(y);
  684.    V_FREE(u);
  685.    V_FREE(v); 
  686.    V_FREE(xn);
  687.    V_FREE(yn);
  688.  
  689.    printf("# Done testing (%s)\n",argv[0]);
  690.    mem_info();
  691. }
  692.