home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / iternsym < prev    next >
Encoding:
Text File  |  1994-02-02  |  29.4 KB  |  1,256 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.c 17/09/93 */
  28.  
  29. /* 
  30.   ITERATIVE METHODS - implementation of several iterative methods;
  31.   see also iter0.c
  32. */
  33.  
  34. #include        <stdio.h>
  35. #include    <math.h>
  36. #include        "matrix.h"
  37. #include        "matrix2.h"
  38. #include    "sparse.h"
  39. #include        "iter.h"
  40.  
  41. static char rcsid[] = "$Id: iternsym.c,v 1.2 1994/02/02 06:02:20 des Exp $";
  42.  
  43.  
  44. #ifdef ANSI_C
  45. VEC    *spCHsolve(SPMAT *,VEC *,VEC *);
  46. #else
  47. VEC    *spCHsolve();
  48. #endif
  49.  
  50.  
  51. /* 
  52.   iter_cgs -- uses CGS to compute a solution x to A.x=b
  53. */
  54.  
  55. VEC    *iter_cgs(ip,r0)
  56. ITER *ip;
  57. VEC *r0;
  58. {
  59.    static VEC  *p = VNULL, *q = VNULL, *r = VNULL, *u = VNULL;
  60.    static VEC  *v = VNULL, *z = VNULL;
  61.    VEC  *tmp;
  62.    Real    alpha, beta, nres, rho, old_rho, sigma, inner;
  63.  
  64.    if (ip == INULL)
  65.      error(E_NULL,"iter_cgs");
  66.    if (!ip->Ax || !ip->b || !r0)
  67.      error(E_NULL,"iter_cgs");
  68.    if ( ip->x == ip->b )
  69.      error(E_INSITU,"iter_cgs");
  70.    if (!ip->stop_crit)
  71.      error(E_NULL,"iter_cgs");
  72.    if ( r0->dim != ip->b->dim )
  73.      error(E_SIZES,"iter_cgs");
  74.    
  75.    if ( ip->eps <= 0.0 )
  76.      ip->eps = MACHEPS;
  77.    
  78.    p = v_resize(p,ip->b->dim);
  79.    q = v_resize(q,ip->b->dim);
  80.    r = v_resize(r,ip->b->dim);
  81.    u = v_resize(u,ip->b->dim);
  82.    v = v_resize(v,ip->b->dim);
  83.  
  84.    MEM_STAT_REG(p,TYPE_VEC);
  85.    MEM_STAT_REG(q,TYPE_VEC);
  86.    MEM_STAT_REG(r,TYPE_VEC);
  87.    MEM_STAT_REG(u,TYPE_VEC);
  88.    MEM_STAT_REG(v,TYPE_VEC);
  89.  
  90.    if (ip->Bx) {
  91.       z = v_resize(z,ip->b->dim);
  92.       MEM_STAT_REG(z,TYPE_VEC); 
  93.    }
  94.  
  95.    if (ip->x != VNULL) {
  96.       if (ip->x->dim != ip->b->dim)
  97.     error(E_SIZES,"iter_cgs");
  98.       ip->Ax(ip->A_par,ip->x,v);            /* v = A*x */
  99.       if (ip->Bx) {
  100.      v_sub(ip->b,v,v);            /* v = b - A*x */
  101.      (ip->Bx)(ip->B_par,v,r);        /* r = B*(b-A*x) */
  102.       }
  103.       else v_sub(ip->b,v,r);            /* r = b-A*x */
  104.    }
  105.    else {  /* ip->x == 0 */
  106.       ip->x = v_get(ip->b->dim);        /* x == 0 */
  107.       ip->shared_x = FALSE;
  108.       if (ip->Bx) (ip->Bx)(ip->B_par,ip->b,r);    /* r = B*b */
  109.       else v_copy(ip->b,r);                       /* r = b */
  110.    }
  111.  
  112.    v_zero(p);    
  113.    v_zero(q);
  114.    old_rho = 1.0;
  115.    
  116.    for (ip->steps = 0; ip->steps <= ip->limit; ip->steps++) {
  117.  
  118.       inner = in_prod(r,r);
  119.       nres = sqrt(fabs(inner));
  120.  
  121.       if (ip->info) ip->info(ip,nres,r,VNULL);
  122.       if ( ip->stop_crit(ip,nres,r,VNULL) ) break;
  123.  
  124.       rho = in_prod(r0,r);
  125.       if ( old_rho == 0.0 )
  126.     error(E_SING,"iter_cgs");
  127.       beta = rho/old_rho;
  128.       v_mltadd(r,q,beta,u);
  129.       v_mltadd(q,p,beta,v);
  130.       v_mltadd(u,v,beta,p);
  131.       
  132.       (ip->Ax)(ip->A_par,p,q);
  133.       if (ip->Bx) {
  134.      (ip->Bx)(ip->B_par,q,z);
  135.      tmp = z;
  136.       }
  137.       else tmp = q;
  138.       
  139.       sigma = in_prod(r0,tmp);
  140.       if ( sigma == 0.0 )
  141.     error(E_SING,"iter_cgs");
  142.       alpha = rho/sigma;
  143.       v_mltadd(u,tmp,-alpha,q);
  144.       v_add(u,q,v);
  145.       
  146.       (ip->Ax)(ip->A_par,v,u);
  147.       if (ip->Bx) {
  148.      (ip->Bx)(ip->B_par,u,z);
  149.      tmp = z;
  150.       }
  151.       else tmp = u;
  152.       
  153.       v_mltadd(r,tmp,-alpha,r);
  154.       v_mltadd(ip->x,v,alpha,ip->x);
  155.       
  156.       old_rho = rho;
  157.    }
  158.  
  159.    return ip->x;
  160. }
  161.  
  162.  
  163.  
  164. /* iter_spcgs -- simple interface for SPMAT data structures 
  165.    use always as follows:
  166.       x = iter_spcgs(A,B,b,r0,tol,x,limit,steps);
  167.    or 
  168.       x = iter_spcgs(A,B,b,r0,tol,VNULL,limit,steps);
  169.    In the second case the solution vector is created.  
  170.    If B is not NULL then it is a preconditioner. 
  171. */
  172. VEC    *iter_spcgs(A,B,b,r0,tol,x,limit,steps)
  173. SPMAT    *A, *B;
  174. VEC    *b, *r0, *x;
  175. double    tol;
  176. int     *steps,limit;
  177. {    
  178.    ITER *ip;
  179.    
  180.    ip = iter_get(0,0);
  181.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  182.    ip->A_par = (void *) A;
  183.    if (B) {
  184.       ip->Bx = (Fun_Ax) sp_mv_mlt;
  185.       ip->B_par = (void *) B;
  186.    }
  187.    else {
  188.       ip->Bx = (Fun_Ax) NULL;
  189.       ip->B_par = NULL;
  190.    }
  191.    ip->info = (Fun_info) NULL;
  192.    ip->limit = limit;
  193.    ip->b = b;
  194.    ip->eps = tol;
  195.    ip->x = x;
  196.    iter_cgs(ip,r0);
  197.    x = ip->x;
  198.    if (steps) *steps = ip->steps;
  199.    ip->shared_x = ip->shared_b = TRUE;   
  200.    iter_free(ip);   /* release only ITER structure */
  201.    return x;        
  202.  
  203. }
  204.  
  205. /*
  206.   Routine for performing LSQR -- the least squares QR algorithm
  207.   of Paige and Saunders:
  208.   "LSQR: an algorithm for sparse linear equations and
  209.   sparse least squares", ACM Trans. Math. Soft., v. 8
  210.   pp. 43--71 (1982)
  211.   */
  212. /* lsqr -- sparse CG-like least squares routine:
  213.    -- finds min_x ||A.x-b||_2 using A defined through A & AT
  214.    -- returns x (if x != NULL) */
  215. VEC    *iter_lsqr(ip)
  216. ITER *ip;
  217. {
  218.    static VEC    *u = VNULL, *v = VNULL, *w = VNULL, *tmp = VNULL;
  219.    Real    alpha, beta, phi, phi_bar;
  220.    Real rho, rho_bar, rho_max, theta, nres;
  221.    Real    s, c;    /* for Givens' rotations */
  222.    int  m, n;
  223.    
  224.    if ( ! ip || ! ip->b || !ip->Ax || !ip->ATx )
  225.      error(E_NULL,"iter_lsqr");
  226.    if ( ip->x == ip->b )
  227.      error(E_INSITU,"iter_lsqr");
  228.    if (!ip->stop_crit || !ip->x)
  229.      error(E_NULL,"iter_lsqr");
  230.  
  231.    if ( ip->eps <= 0.0 )
  232.      ip->eps = MACHEPS;
  233.    
  234.    m = ip->b->dim;    
  235.    n = ip->x->dim;
  236.  
  237.    u = v_resize(u,(u_int)m);
  238.    v = v_resize(v,(u_int)n);
  239.    w = v_resize(w,(u_int)n);
  240.    tmp = v_resize(tmp,(u_int)n);
  241.  
  242.    MEM_STAT_REG(u,TYPE_VEC);
  243.    MEM_STAT_REG(v,TYPE_VEC);
  244.    MEM_STAT_REG(w,TYPE_VEC);
  245.    MEM_STAT_REG(tmp,TYPE_VEC);  
  246.  
  247.    if (ip->x != VNULL) {
  248.       ip->Ax(ip->A_par,ip->x,u);            /* u = A*x */
  249.       v_sub(ip->b,u,u);                /* u = b-A*x */
  250.    }
  251.    else {  /* ip->x == 0 */
  252.       ip->x = v_get(ip->b->dim);
  253.       ip->shared_x = FALSE;
  254.       v_copy(ip->b,u);                       /* u = b */
  255.    }
  256.  
  257.    beta = v_norm2(u); 
  258.    if ( beta == 0.0 )
  259.      return ip->x;
  260.    sv_mlt(1.0/beta,u,u);
  261.    (ip->ATx)(ip->AT_par,u,v);
  262.    alpha = v_norm2(v);
  263.    if ( alpha == 0.0 )
  264.      return ip->x;
  265.    sv_mlt(1.0/alpha,v,v);
  266.    v_copy(v,w);
  267.    phi_bar = beta;
  268.    rho_bar = alpha;
  269.    
  270.    rho_max = 1.0;
  271.    for (ip->steps = 0; ip->steps <= ip->limit; ip->steps++) {
  272.  
  273.       tmp = v_resize(tmp,m);
  274.       (ip->Ax)(ip->A_par,v,tmp);
  275.       
  276.       v_mltadd(tmp,u,-alpha,u);
  277.       beta = v_norm2(u);    
  278.       sv_mlt(1.0/beta,u,u);
  279.       
  280.       tmp = v_resize(tmp,n);
  281.       (ip->ATx)(ip->AT_par,u,tmp);
  282.       v_mltadd(tmp,v,-beta,v);
  283.       alpha = v_norm2(v);    
  284.       sv_mlt(1.0/alpha,v,v);
  285.       
  286.       rho = sqrt(rho_bar*rho_bar+beta*beta);
  287.       if ( rho > rho_max )
  288.     rho_max = rho;
  289.       c   = rho_bar/rho;
  290.       s   = beta/rho;
  291.       theta   =  s*alpha;
  292.       rho_bar = -c*alpha;
  293.       phi     =  c*phi_bar;
  294.       phi_bar =  s*phi_bar;
  295.       
  296.       /* update ip->x & w */
  297.       if ( rho == 0.0 )
  298.     error(E_SING,"iter_lsqr");
  299.       v_mltadd(ip->x,w,phi/rho,ip->x);
  300.       v_mltadd(v,w,-theta/rho,w);
  301.  
  302.       nres = fabs(phi_bar*alpha*c)*rho_max;
  303.  
  304.       if (ip->info) ip->info(ip,nres,w,VNULL);
  305.       if ( ip->stop_crit(ip,nres,w,VNULL) ) break;
  306.    } 
  307.    
  308.    return ip->x;
  309. }
  310.  
  311. /* iter_splsqr -- simple interface for SPMAT data structures */
  312. VEC    *iter_splsqr(A,b,tol,x,limit,steps)
  313. SPMAT    *A;
  314. VEC    *b, *x;
  315. double    tol;
  316. int *steps,limit;
  317. {
  318.    ITER *ip;
  319.    
  320.    ip = iter_get(0,0);
  321.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  322.    ip->A_par = (void *) A;
  323.    ip->ATx = (Fun_Ax) sp_vm_mlt;
  324.    ip->AT_par = (void *) A;
  325.    ip->Bx = (Fun_Ax) NULL;
  326.    ip->B_par = NULL;
  327.  
  328.    ip->info = (Fun_info) NULL;
  329.    ip->limit = limit;
  330.    ip->b = b;
  331.    ip->eps = tol;
  332.    ip->x = x;
  333.    iter_lsqr(ip);
  334.    x = ip->x;
  335.    if (steps) *steps = ip->steps;
  336.    ip->shared_x = ip->shared_b = TRUE;
  337.    iter_free(ip);   /* release only ITER structure */
  338.    return x;        
  339. }
  340.  
  341.  
  342.  
  343. /* iter_arnoldi -- an implementation of the Arnoldi method;
  344.    iterative refinement is applied.
  345. */
  346. MAT    *iter_arnoldi_iref(ip,h_rem,Q,H)
  347. ITER  *ip;
  348. Real  *h_rem;
  349. MAT   *Q, *H;
  350. {
  351.    static VEC *u=VNULL, *r=VNULL, *s=VNULL, *tmp=VNULL;
  352.    VEC v;     /* auxiliary vector */
  353.    int    i,j;
  354.    Real    h_val, c;
  355.    
  356.    if (ip == INULL)
  357.      error(E_NULL,"iter_arnoldi_iref");
  358.    if ( ! ip->Ax || ! Q || ! ip->x )
  359.      error(E_NULL,"iter_arnoldi_iref");
  360.    if ( ip->k <= 0 )
  361.      error(E_BOUNDS,"iter_arnoldi_iref");
  362.    if ( Q->n != ip->x->dim ||    Q->m != ip->k )
  363.      error(E_SIZES,"iter_arnoldi_iref");
  364.    
  365.    m_zero(Q);
  366.    H = m_resize(H,ip->k,ip->k);
  367.    m_zero(H);
  368.  
  369.    u = v_resize(u,ip->x->dim);
  370.    tmp = v_resize(tmp,ip->x->dim);
  371.    r = v_resize(r,ip->k);
  372.    s = v_resize(s,ip->k);
  373.    MEM_STAT_REG(u,TYPE_VEC);
  374.    MEM_STAT_REG(tmp,TYPE_VEC);
  375.    MEM_STAT_REG(r,TYPE_VEC);
  376.    MEM_STAT_REG(s,TYPE_VEC);
  377.  
  378.    v.dim = v.max_dim = ip->x->dim;
  379.  
  380.    c = v_norm2(ip->x);
  381.    if ( c <= 0.0)
  382.      return H;
  383.    else {
  384.       v.ve = Q->me[0];
  385.       sv_mlt(1.0/c,ip->x,&v);
  386.    }
  387.  
  388.    v_zero(r);
  389.    v_zero(s);
  390.    for ( i = 0; i < ip->k; i++ )
  391.    {
  392.       v.ve = Q->me[i];
  393.       u = (ip->Ax)(ip->A_par,&v,u);
  394.       v_zero(tmp);
  395.       for (j = 0; j <= i; j++) {
  396.      v.ve = Q->me[j];
  397.      r->ve[j] = in_prod(&v,u);
  398.      v_mltadd(tmp,&v,r->ve[j],tmp);
  399.       }
  400.       v_sub(u,tmp,u);
  401.       h_val = v_norm2(u);
  402.       /* if u == 0 then we have an exact subspace */
  403.       if ( h_val <= 0.0 )
  404.       {
  405.      *h_rem = h_val;
  406.      return H;
  407.       }
  408.       /* iterative refinement -- ensures near orthogonality */
  409.       do {
  410.      v_zero(tmp);
  411.      for (j = 0; j <= i; j++) {
  412.         v.ve = Q->me[j];
  413.         s->ve[j] = in_prod(&v,u);
  414.         v_mltadd(tmp,&v,s->ve[j],tmp);
  415.      }
  416.      v_sub(u,tmp,u);
  417.      v_add(r,s,r);
  418.       } while ( v_norm2(s) > 0.1*(h_val = v_norm2(u)) );
  419.       /* now that u is nearly orthogonal to Q, update H */
  420.       set_col(H,i,r);
  421.       /* check once again if h_val is zero */
  422.       if ( h_val <= 0.0 )
  423.       {
  424.      *h_rem = h_val;
  425.      return H;
  426.       }
  427.       if ( i == ip->k-1 )
  428.       {
  429.      *h_rem = h_val;
  430.      continue;
  431.       }
  432.       /* H->me[i+1][i] = h_val; */
  433.       m_set_val(H,i+1,i,h_val);
  434.       v.ve = Q->me[i+1];
  435.       sv_mlt(1.0/h_val,u,&v);
  436.    }
  437.    
  438.    return H;
  439. }
  440.  
  441. /* iter_arnoldi -- an implementation of the Arnoldi method;
  442.    without iterative refinement
  443. */
  444. MAT    *iter_arnoldi(ip,h_rem,Q,H)
  445. ITER  *ip;
  446. Real  *h_rem;
  447. MAT   *Q, *H;
  448. {
  449.    static VEC *u=VNULL, *r=VNULL, *tmp=VNULL;
  450.    VEC v;     /* auxiliary vector */
  451.    int    i,j;
  452.    Real    h_val, c;
  453.    
  454.    if (ip == INULL)
  455.      error(E_NULL,"iter_arnoldi");
  456.    if ( ! ip->Ax || ! Q || ! ip->x )
  457.      error(E_NULL,"iter_arnoldi");
  458.    if ( ip->k <= 0 )
  459.      error(E_BOUNDS,"iter_arnoldi");
  460.    if ( Q->n != ip->x->dim ||    Q->m != ip->k )
  461.      error(E_SIZES,"iter_arnoldi");
  462.    
  463.    m_zero(Q);
  464.    H = m_resize(H,ip->k,ip->k);
  465.    m_zero(H);
  466.  
  467.    u = v_resize(u,ip->x->dim);
  468.    tmp = v_resize(tmp,ip->x->dim);
  469.    r = v_resize(r,ip->k);
  470.    MEM_STAT_REG(u,TYPE_VEC);
  471.    MEM_STAT_REG(tmp,TYPE_VEC);
  472.    MEM_STAT_REG(r,TYPE_VEC);
  473.  
  474.    v.dim = v.max_dim = ip->x->dim;
  475.  
  476.    c = v_norm2(ip->x);
  477.    if ( c <= 0.0)
  478.      return H;
  479.    else {
  480.       v.ve = Q->me[0];
  481.       sv_mlt(1.0/c,ip->x,&v);
  482.    }
  483.  
  484.    v_zero(r);
  485.    for ( i = 0; i < ip->k; i++ )
  486.    {
  487.       v.ve = Q->me[i];
  488.       u = (ip->Ax)(ip->A_par,&v,u);
  489.       v_zero(tmp);
  490.       for (j = 0; j <= i; j++) {
  491.      v.ve = Q->me[j];
  492.      r->ve[j] = in_prod(&v,u);
  493.      v_mltadd(tmp,&v,r->ve[j],tmp);
  494.       }
  495.       v_sub(u,tmp,u);
  496.       h_val = v_norm2(u);
  497.       /* if u == 0 then we have an exact subspace */
  498.       if ( h_val <= 0.0 )
  499.       {
  500.      *h_rem = h_val;
  501.      return H;
  502.       }
  503.       set_col(H,i,r);
  504.       if ( i == ip->k-1 )
  505.       {
  506.      *h_rem = h_val;
  507.      continue;
  508.       }
  509.       /* H->me[i+1][i] = h_val; */
  510.       m_set_val(H,i+1,i,h_val);
  511.       v.ve = Q->me[i+1];
  512.       sv_mlt(1.0/h_val,u,&v);
  513.    }
  514.    
  515.    return H;
  516. }
  517.  
  518.  
  519.  
  520. /* iter_sparnoldi -- uses arnoldi() with an explicit representation of A */
  521. MAT    *iter_sparnoldi(A,x0,m,h_rem,Q,H)
  522. SPMAT    *A;
  523. VEC    *x0;
  524. int    m;
  525. Real    *h_rem;
  526. MAT    *Q, *H;
  527. {
  528.    ITER *ip;
  529.    
  530.    ip = iter_get(0,0);
  531.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  532.    ip->A_par = (void *) A;
  533.    ip->x = x0;
  534.    ip->k = m;
  535.    iter_arnoldi_iref(ip,h_rem,Q,H);
  536.    ip->shared_x = ip->shared_b = TRUE;
  537.    iter_free(ip);   /* release only ITER structure */
  538.    return H;    
  539. }
  540.  
  541.  
  542. /* for testing gmres */
  543. static void test_gmres(ip,i,Q,R,givc,givs,h_val)
  544. ITER *ip;
  545. int i;
  546. MAT *Q, *R;
  547. VEC *givc, *givs;
  548. double h_val;
  549. {
  550.    VEC vt, vt1;
  551.    static MAT *Q1, *R1;
  552.    int j;
  553.    
  554.    /* test Q*A*Q^T = R  */
  555.  
  556.    Q = m_resize(Q,i+1,ip->b->dim);
  557.    Q1 = m_resize(Q1,i+1,ip->b->dim);
  558.    R1 = m_resize(R1,i+1,i+1);
  559.    MEM_STAT_REG(Q1,TYPE_MAT);
  560.    MEM_STAT_REG(R1,TYPE_MAT);
  561.  
  562.    vt.dim = vt.max_dim = ip->b->dim;
  563.    vt1.dim = vt1.max_dim = ip->b->dim;
  564.    for (j=0; j <= i; j++) {
  565.       vt.ve = Q->me[j];
  566.       vt1.ve = Q1->me[j];
  567.       ip->Ax(ip->A_par,&vt,&vt1);
  568.    }
  569.  
  570.    mmtr_mlt(Q,Q1,R1);
  571.    R1 = m_resize(R1,i+2,i+1);
  572.    for (j=0; j < i; j++)
  573.      R1->me[i+1][j] = 0.0;
  574.    R1->me[i+1][i] = h_val;
  575.    
  576.    for (j = 0; j <= i; j++) {
  577.       rot_rows(R1,j,j+1,givc->ve[j],givs->ve[j],R1);
  578.    }
  579.  
  580.    R1 = m_resize(R1,i+1,i+1);
  581.    m_sub(R,R1,R1);
  582.    /* if (m_norm_inf(R1) > MACHEPS*ip->b->dim)  */
  583.    printf(" %d. ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
  584.       ip->steps,m_norm_inf(R1),MACHEPS);
  585.    
  586.    /* check Q*Q^T = I */
  587.    
  588.    Q = m_resize(Q,i+1,ip->b->dim);
  589.    mmtr_mlt(Q,Q,R1);
  590.    for (j=0; j <= i; j++)
  591.      R1->me[j][j] -= 1.0;
  592.    if (m_norm_inf(R1) > MACHEPS*ip->b->dim)
  593.      printf(" ! m_norm_inf(Q*Q^T) = %g\n",m_norm_inf(R1));  
  594.    
  595. }
  596.  
  597.  
  598. /* gmres -- generalised minimum residual algorithm of Saad & Schultz
  599.    SIAM J. Sci. Stat. Comp. v.7, pp.856--869 (1986)
  600. */
  601. VEC    *iter_gmres(ip)
  602. ITER *ip;
  603. {
  604.    static VEC *u=VNULL, *r=VNULL, *rhs = VNULL;
  605.    static VEC *givs=VNULL, *givc=VNULL, *z = VNULL;
  606.    static MAT *Q = MNULL, *R = MNULL;
  607.    VEC *rr, v, v1;   /* additional pointers (not real vectors) */
  608.    int    i,j, done;
  609.    Real    nres;
  610. /*   Real last_h;  */
  611.    
  612.    if (ip == INULL)
  613.      error(E_NULL,"iter_gmres");
  614.    if ( ! ip->Ax || ! ip->b )
  615.      error(E_NULL,"iter_gmres");
  616.    if ( ! ip->stop_crit )
  617.      error(E_NULL,"iter_gmres");
  618.    if ( ip->k <= 0 )
  619.      error(E_BOUNDS,"iter_gmres");
  620.    if (ip->x != VNULL && ip->x->dim != ip->b->dim)
  621.      error(E_SIZES,"iter_gmres");
  622.  
  623.    r = v_resize(r,ip->k+1);
  624.    u = v_resize(u,ip->b->dim);
  625.    rhs = v_resize(rhs,ip->k+1);
  626.    givs = v_resize(givs,ip->k);  /* Givens rotations */
  627.    givc = v_resize(givc,ip->k); 
  628.    
  629.    MEM_STAT_REG(r,TYPE_VEC);
  630.    MEM_STAT_REG(u,TYPE_VEC);
  631.    MEM_STAT_REG(rhs,TYPE_VEC);
  632.    MEM_STAT_REG(givs,TYPE_VEC);
  633.    MEM_STAT_REG(givc,TYPE_VEC);
  634.    
  635.    R = m_resize(R,ip->k+1,ip->k);
  636.    Q = m_resize(Q,ip->k,ip->b->dim);
  637.    MEM_STAT_REG(R,TYPE_MAT);
  638.    MEM_STAT_REG(Q,TYPE_MAT);        
  639.  
  640.    if (ip->x == VNULL) {  /* ip->x == 0 */
  641.       ip->x = v_get(ip->b->dim);
  642.       ip->shared_x = FALSE;
  643.    }   
  644.  
  645.    v.dim = v.max_dim = ip->b->dim;      /* v and v1 are pointers to rows */
  646.    v1.dim = v1.max_dim = ip->b->dim;      /* of matrix Q */
  647.    
  648.    if (ip->Bx != (Fun_Ax)NULL) {    /* if precondition is defined */
  649.       z = v_resize(z,ip->b->dim);
  650.       MEM_STAT_REG(z,TYPE_VEC);
  651.    }
  652.    
  653.    done = FALSE;
  654.    for (ip->steps = 0; ip->steps <= ip->limit; ) {
  655.  
  656.       /* restart */
  657.  
  658.       ip->Ax(ip->A_par,ip->x,u);            /* u = A*x */
  659.       v_sub(ip->b,u,u);                 /* u = b - A*x */
  660.       rr = u;                /* rr is a pointer only */
  661.       
  662.       if (ip->Bx) {
  663.      (ip->Bx)(ip->B_par,u,z);            /* tmp = B*(b-A*x)  */
  664.      rr = z;
  665.       }
  666.       
  667.       nres = v_norm2(rr);
  668.       if (ip->steps == 0) {
  669.      if (ip->info) ip->info(ip,nres,VNULL,VNULL);
  670.      if ( ip->stop_crit(ip,nres,VNULL,VNULL) ) {
  671.         done = TRUE;
  672.         break;
  673.      }
  674.       }
  675.       else if (nres <= 0.0) break;
  676.  
  677.       v.ve = Q->me[0];
  678.       sv_mlt(1.0/nres,rr,&v);
  679.       
  680.       v_zero(r);
  681.       v_zero(rhs);
  682.       rhs->ve[0] = nres;
  683.  
  684.       for ( i = 0; i < ip->k; i++ ) {
  685.      ip->steps++;
  686.      v.ve = Q->me[i];    
  687.      (ip->Ax)(ip->A_par,&v,u);
  688.      rr = u;
  689.      if (ip->Bx) {
  690.         (ip->Bx)(ip->B_par,u,z);
  691.         rr = z;
  692.      }
  693.      
  694.      if (i < ip->k - 1) {
  695.         v1.ve = Q->me[i+1];
  696.         v_copy(rr,&v1);
  697.         for (j = 0; j <= i; j++) {
  698.            v.ve = Q->me[j];
  699.            r->ve[j] = in_prod(&v,rr);
  700.            v_mltadd(&v1,&v,-r->ve[j],&v1);
  701.         }
  702.         
  703.         r->ve[i+1] = nres = v_norm2(&v1);
  704.         if (nres <= 0.0) {
  705.            warning(WARN_RES_LESS_0,"iter_gmres");
  706.            break;
  707.         }
  708.         sv_mlt(1.0/nres,&v1,&v1);
  709.      }
  710.      else {  /* i == ip->k - 1 */
  711.         /* Q->me[ip->k] need not be computed */
  712.  
  713.         for (j = 0; j <= i; j++) {
  714.            v.ve = Q->me[j];
  715.            r->ve[j] = in_prod(&v,rr);
  716.         }
  717.         
  718.         nres = in_prod(rr,rr) - in_prod(r,r);
  719.         if (nres <= 0.0) {
  720.            warning(WARN_RES_LESS_0,"iter_gmres");
  721.            break;
  722.         }
  723.         r->ve[i+1] = sqrt(nres);
  724.      }
  725.  
  726.      /* QR update */
  727.  
  728.      /* last_h = r->ve[i+1]; */ /* for test only */
  729.      for (j = 0; j < i; j++) 
  730.        rot_vec(r,j,j+1,givc->ve[j],givs->ve[j],r);
  731.      givens(r->ve[i],r->ve[i+1],&givc->ve[i],&givs->ve[i]);
  732.      rot_vec(r,i,i+1,givc->ve[i],givs->ve[i],r);
  733.      rot_vec(rhs,i,i+1,givc->ve[i],givs->ve[i],rhs);
  734.      
  735.      set_col(R,i,r);
  736.  
  737.      nres = fabs((double) rhs->ve[i+1]);
  738.      if (ip->info) ip->info(ip,nres,VNULL,VNULL);
  739.      if ( ip->stop_crit(ip,nres,VNULL,VNULL) ||
  740.          ip->steps >= ip->limit ) {
  741.         done = TRUE;
  742.         break;
  743.      }
  744.       }
  745.       
  746.       /* use ixi submatrix of R */
  747.  
  748.       if (nres <= 0.0) {
  749.      i--;
  750.      done = TRUE;
  751.       }
  752.       if (i == ip->k) i = ip->k - 1;
  753.  
  754.       R = m_resize(R,i+1,i+1);
  755.       rhs = v_resize(rhs,i+1);
  756.       
  757.       /* test only */
  758.       /* test_gmres(ip,i,Q,R,givc,givs,last_h);  */
  759.       
  760.       Usolve(R,rhs,rhs,0.0);      /* solve a system: R*x = rhs */
  761.  
  762.       /* new approximation */
  763.  
  764.       for (j = 0; j <= i; j++) {
  765.      v.ve = Q->me[j]; 
  766.      v_mltadd(ip->x,&v,rhs->ve[j],ip->x);
  767.       }
  768.  
  769.       if (done) break;
  770.  
  771.       /* back to old dimensions */
  772.  
  773.       rhs = v_resize(rhs,ip->k+1);
  774.       R = m_resize(R,ip->k+1,ip->k);
  775.  
  776.    }
  777.  
  778.    return ip->x;
  779. }
  780.  
  781. /* iter_spgmres - a simple interface to iter_gmres */
  782.  
  783. VEC    *iter_spgmres(A,B,b,tol,x,k,limit,steps)
  784. SPMAT    *A, *B;
  785. VEC    *b, *x;
  786. double    tol;
  787. int *steps,k,limit;
  788. {
  789.    ITER *ip;
  790.    
  791.    ip = iter_get(0,0);
  792.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  793.    ip->A_par = (void *) A;
  794.    if (B) {
  795.       ip->Bx = (Fun_Ax) sp_mv_mlt;
  796.       ip->B_par = (void *) B;
  797.    }
  798.    else {
  799.       ip->Bx = (Fun_Ax) NULL;
  800.       ip->B_par = NULL;
  801.    }
  802.    ip->k = k;
  803.    ip->limit = limit;
  804.    ip->info = (Fun_info) NULL;
  805.    ip->b = b;
  806.    ip->eps = tol;
  807.    ip->x = x;
  808.    iter_gmres(ip);
  809.    x = ip->x;
  810.    if (steps) *steps = ip->steps;
  811.    ip->shared_x = ip->shared_b = TRUE;
  812.    iter_free(ip);   /* release only ITER structure */
  813.    return x;        
  814. }
  815.  
  816.  
  817. /* for testing mgcr */
  818. static void test_mgcr(ip,i,Q,R)
  819. ITER *ip;
  820. int i;
  821. MAT *Q, *R;
  822. {
  823.    VEC vt, vt1;
  824.    static MAT *R1;
  825.    static VEC *r, *r1;
  826.    VEC *rr;
  827.    int k,j;
  828.    Real sm;
  829.    
  830.    
  831.    /* check Q*Q^T = I */
  832.    vt.dim = vt.max_dim = ip->b->dim;
  833.    vt1.dim = vt1.max_dim = ip->b->dim;
  834.    
  835.    Q = m_resize(Q,i+1,ip->b->dim);
  836.    R1 = m_resize(R1,i+1,i+1);
  837.    r = v_resize(r,ip->b->dim);
  838.    r1 = v_resize(r1,ip->b->dim);
  839.    MEM_STAT_REG(R1,TYPE_MAT);
  840.    MEM_STAT_REG(r,TYPE_VEC);
  841.    MEM_STAT_REG(r1,TYPE_VEC);
  842.  
  843.    m_zero(R1);
  844.    for (k=1; k <= i; k++)
  845.      for (j=1; j <= i; j++) {
  846.     vt.ve = Q->me[k];
  847.     vt1.ve = Q->me[j];
  848.     R1->me[k][j] = in_prod(&vt,&vt1);
  849.      }
  850.    for (j=1; j <= i; j++)
  851.      R1->me[j][j] -= 1.0;
  852.    if (m_norm_inf(R1) > MACHEPS*ip->b->dim)
  853.      printf(" ! (mgcr:) m_norm_inf(Q*Q^T) = %g\n",m_norm_inf(R1));  
  854.  
  855.    /* check (r_i,Ap_j) = 0 for j <= i */
  856.    
  857.    ip->Ax(ip->A_par,ip->x,r);
  858.    v_sub(ip->b,r,r);
  859.    rr = r;
  860.    if (ip->Bx) {
  861.       ip->Bx(ip->B_par,r,r1);
  862.       rr = r1;
  863.    }
  864.    
  865.    printf(" ||r|| = %g\n",v_norm2(rr));
  866.    sm = 0.0;
  867.    for (j = 1; j <= i; j++) {
  868.       vt.ve = Q->me[j];
  869.       sm = max(sm,in_prod(&vt,rr));
  870.    }
  871.    if (sm >= MACHEPS*ip->b->dim)
  872.      printf(" ! (mgcr:) max_j (r,Ap_j) = %g\n",sm);
  873.  
  874. }
  875.  
  876.  
  877.  
  878.  
  879. /* 
  880.   iter_mgcr -- modified generalized conjugate residual algorithm;
  881.   fast version of GCR;
  882. */
  883. VEC *iter_mgcr(ip)
  884. ITER *ip;
  885. {
  886.    static VEC *As, *beta, *alpha, *z;
  887.    static MAT *N, *H;
  888.    
  889.    VEC *rr, v, s;  /* additional pointer and structures */
  890.    Real nres;      /* norm of a residual */
  891.    Real dd;        /* coefficient d_i */
  892.    int i,j;
  893.    int done;      /* if TRUE then stop the iterative process */
  894.    int dim;       /* dimension of the problem */
  895.    
  896.    /* ip cannot be NULL */
  897.    if (ip == INULL) error(E_NULL,"mgcr");
  898.    /* Ax, b and stopping criterion must be given */
  899.    if (! ip->Ax || ! ip->b || ! ip->stop_crit) 
  900.      error(E_NULL,"mgcr");
  901.    /* at least one direction vector must exist */
  902.    if ( ip->k <= 0) error(E_BOUNDS,"mgcr");
  903.    /* if the vector x is given then b and x must have the same dimension */
  904.    if ( ip->x && ip->x->dim != ip->b->dim)
  905.      error(E_SIZES,"mgcr");
  906.    
  907.    dim = ip->b->dim;
  908.    As = v_resize(As,dim);
  909.    alpha = v_resize(alpha,ip->k);
  910.    beta = v_resize(beta,ip->k);
  911.    
  912.    MEM_STAT_REG(As,TYPE_VEC);
  913.    MEM_STAT_REG(alpha,TYPE_VEC);
  914.    MEM_STAT_REG(beta,TYPE_VEC);
  915.    
  916.    H = m_resize(H,ip->k,ip->k);
  917.    N = m_resize(N,ip->k,dim);
  918.    
  919.    MEM_STAT_REG(H,TYPE_MAT);
  920.    MEM_STAT_REG(N,TYPE_MAT);
  921.    
  922.    /* if a preconditioner is defined */
  923.    if (ip->Bx) {
  924.       z = v_resize(z,dim);
  925.       MEM_STAT_REG(z,TYPE_VEC);
  926.    }
  927.    
  928.    /* if x is NULL then it is assumed that x has 
  929.       entries with value zero */
  930.    if ( ! ip->x ) {
  931.       ip->x = v_get(ip->b->dim);
  932.       ip->shared_x = FALSE;
  933.    }
  934.    
  935.    /* v and s are additional pointers to rows of N */
  936.    /* they must have the same dimension as rows of N */
  937.    v.dim = v.max_dim = s.dim = s.max_dim = dim;
  938.    
  939.    
  940.    done = FALSE;
  941.    ip->steps = 0;
  942.    while (TRUE) {
  943.       (*ip->Ax)(ip->A_par,ip->x,As);         /* As = A*x */
  944.       v_sub(ip->b,As,As);                    /* As = b - A*x */
  945.       rr = As;                               /* rr is an additional pointer */
  946.       
  947.       /* if a preconditioner is defined */
  948.       if (ip->Bx) {
  949.      (*ip->Bx)(ip->B_par,As,z);               /* z = B*(b-A*x)  */
  950.      rr = z;                                  
  951.       }
  952.       
  953.       /* norm of the residual */
  954.       nres = v_norm2(rr);
  955.       dd = nres*nres;                            /* dd = ||r_i||^2  */
  956.       
  957.       /* we need to check if the norm of the residual is small enough 
  958.      only when we start the iterative process;
  959.      otherwise it has been checked yet. 
  960.      Also the member ip->init_res is updated indirectly by 
  961.      ip->stop_crit.
  962.      */
  963.       if (ip->steps == 0) {                /* information for a user */
  964.      if (ip->info) (*ip->info)(ip,nres,As,rr); 
  965.      if ( (*ip->stop_crit)(ip,nres,As,rr) ) { 
  966.           /* iterative process is finished */
  967.         done = TRUE; 
  968.         break;
  969.      }
  970.       }
  971.       else if (nres <= 0.0) 
  972.     break;  /* residual is zero -> finish the process */ 
  973.       
  974.       /* save this residual in the first row of N */
  975.       v.ve = N->me[0];
  976.       v_copy(rr,&v);
  977.       
  978.       for (i = 0; i < ip->k && ip->steps <= ip->limit; i++) {
  979.      ip->steps++;
  980.      v.ve = N->me[i];                /* pointer to a row of N (=s_i) */
  981.      /* note that we must use here &v, not v */
  982.      (*ip->Ax)(ip->A_par,&v,As); 
  983.      rr = As;                        /* As = A*s_i */
  984.      if (ip->Bx) {
  985.         (*ip->Bx)(ip->B_par,As,z);    /* z = B*A*s_i  */
  986.         rr = z;
  987.      }
  988.      
  989.      if (i < ip->k - 1) {
  990.         s.ve = N->me[i+1];         /* pointer to a row of N (=s_{i+1}) */
  991.         v_copy(rr,&s);                   /* s_{i+1} = B*A*s_i */
  992.         for (j = 0; j <= i-1; j++) {
  993.            v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
  994.            beta->ve[j] = in_prod(&v,rr);    /* beta_{j,i} */
  995.                    /* s_{i+1} -= beta_{j,i}*s_{j+1} */
  996.            v_mltadd(&s,&v,- beta->ve[j],&s);    
  997.         }
  998.         
  999.          /* beta_{i,i} = ||s_{i+1}||_2 */
  1000.         beta->ve[i] = nres = v_norm2(&s);     
  1001.         if ( nres <= 0.0) break;         /* if s_{i+1} == 0 then finish */
  1002.         sv_mlt(1.0/nres,&s,&s);           /* normalize s_{i+1} */
  1003.         
  1004.         v.ve = N->me[0];
  1005.         alpha->ve[i] = in_prod(&v,&s);     /* alpha_i = (s_0 , s_{i+1}) */
  1006.         
  1007.      }
  1008.      else {
  1009.         for (j = 0; j <= i-1; j++) {
  1010.            v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
  1011.            beta->ve[j] = in_prod(&v,rr);       /* beta_{j,i} */
  1012.         }
  1013.         
  1014.         nres = in_prod(rr,rr);                 /* rr = B*A*s_{k-1} */
  1015.         for (j = 0; j <= i-1; j++)
  1016.               nres -= beta->ve[j]*beta->ve[j];
  1017.         if (nres <= 0.0)  break;               /* s_k is zero */
  1018.         else beta->ve[i] = sqrt(nres);         /* beta_{k-1,k-1} */
  1019.         
  1020.         v.ve = N->me[0];
  1021.         alpha->ve[i] = in_prod(&v,rr); 
  1022.         for (j = 0; j <= i-1; j++)
  1023.               alpha->ve[i] -= beta->ve[j]*alpha->ve[j];
  1024.         alpha->ve[i] /= beta->ve[i];                /* alpha_{k-1} */
  1025.         
  1026.      }
  1027.      
  1028.      set_col(H,i,beta);
  1029.      
  1030.      dd -= alpha->ve[i]*alpha->ve[i];
  1031.      nres = sqrt(fabs((double) dd));
  1032.      if (dd < 0.0)  {     /* do restart */
  1033.         if (ip->info) (*ip->info)(ip,-nres,VNULL,VNULL);  
  1034.         break;
  1035.      }
  1036.      
  1037.      if (ip->info) (*ip->info)(ip,nres,VNULL,VNULL);     
  1038.      if ( ip->stop_crit(ip,nres,VNULL,VNULL) ) {
  1039.         /* stopping criterion is satisfied */
  1040.         done = TRUE;
  1041.         break;
  1042.      }
  1043.      
  1044.       } /* end of for */
  1045.       
  1046.       if (nres <= 0.0) {
  1047.      i--;
  1048.      done = TRUE;
  1049.       }
  1050.       if (i >= ip->k) i = ip->k - 1;
  1051.       
  1052.       /* use (i+1) by (i+1) submatrix of H */
  1053.       H = m_resize(H,i+1,i+1);
  1054.       alpha = v_resize(alpha,i+1);
  1055.       Usolve(H,alpha,alpha,0.0);       /* c_i is saved in alpha */
  1056.       
  1057.       for (j = 0; j <= i; j++) {
  1058.      v.ve = N->me[j];
  1059.      v_mltadd(ip->x,&v,alpha->ve[j],ip->x);
  1060.       }
  1061.       
  1062.       
  1063.       if (done) break;              /* stop the iterative process */
  1064.       alpha = v_resize(alpha,ip->k);
  1065.       H = m_resize(H,ip->k,ip->k);
  1066.       
  1067.    }  /* end of while */
  1068.    
  1069.    return ip->x;                    /* return the solution */
  1070. }
  1071.  
  1072.  
  1073.  
  1074. /* iter_spmgcr - a simple interface to iter_mgcr */
  1075. /* no preconditioner */
  1076. VEC    *iter_spmgcr(A,B,b,tol,x,k,limit,steps)
  1077. SPMAT    *A, *B;
  1078. VEC    *b, *x;
  1079. double    tol;
  1080. int *steps,k,limit;
  1081. {
  1082.    ITER *ip;
  1083.    
  1084.    ip = iter_get(0,0);
  1085.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  1086.    ip->A_par = (void *) A;
  1087.    if (B) {
  1088.       ip->Bx = (Fun_Ax) sp_mv_mlt;
  1089.       ip->B_par = (void *) B;
  1090.    }
  1091.    else {
  1092.       ip->Bx = (Fun_Ax) NULL;
  1093.       ip->B_par = NULL;
  1094.    }
  1095.  
  1096.    ip->k = k;
  1097.    ip->limit = limit;
  1098.    ip->info = (Fun_info) NULL;
  1099.    ip->b = b;
  1100.    ip->eps = tol;
  1101.    ip->x = x;
  1102.    iter_mgcr(ip);
  1103.    x = ip->x;
  1104.    if (steps) *steps = ip->steps;
  1105.    ip->shared_x = ip->shared_b = TRUE;
  1106.    iter_free(ip);   /* release only ITER structure */
  1107.    return x;        
  1108. }
  1109.  
  1110.  
  1111.  
  1112. /* 
  1113.   Conjugate gradients method for a normal equation
  1114.   a preconditioner B must be symmetric !!
  1115. */
  1116. VEC  *iter_cgne(ip)
  1117. ITER *ip;
  1118. {
  1119.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  1120.    Real    alpha, beta, inner, old_inner, nres;
  1121.    VEC *rr1;   /* pointer only */
  1122.    
  1123.    if (ip == INULL)
  1124.      error(E_NULL,"iter_cgne");
  1125.    if (!ip->Ax || ! ip->ATx || !ip->b)
  1126.      error(E_NULL,"iter_cgne");
  1127.    if ( ip->x == ip->b )
  1128.      error(E_INSITU,"iter_cgne");
  1129.    if (!ip->stop_crit)
  1130.      error(E_NULL,"iter_cgne");
  1131.    
  1132.    if ( ip->eps <= 0.0 )
  1133.      ip->eps = MACHEPS;
  1134.    
  1135.    r = v_resize(r,ip->b->dim);
  1136.    p = v_resize(p,ip->b->dim);
  1137.    q = v_resize(q,ip->b->dim);
  1138.  
  1139.    MEM_STAT_REG(r,TYPE_VEC);
  1140.    MEM_STAT_REG(p,TYPE_VEC);
  1141.    MEM_STAT_REG(q,TYPE_VEC);
  1142.  
  1143.    z = v_resize(z,ip->b->dim);
  1144.    MEM_STAT_REG(z,TYPE_VEC);
  1145.  
  1146.    if (ip->x) {
  1147.       if (ip->x->dim != ip->b->dim)
  1148.     error(E_SIZES,"iter_cgne");
  1149.       ip->Ax(ip->A_par,ip->x,p);            /* p = A*x */
  1150.       v_sub(ip->b,p,z);                 /* z = b - A*x */
  1151.    }
  1152.    else {  /* ip->x == 0 */
  1153.       ip->x = v_get(ip->b->dim);
  1154.       ip->shared_x = FALSE;
  1155.       v_copy(ip->b,z);
  1156.    }
  1157.    rr1 = z;
  1158.    if (ip->Bx) {
  1159.       (ip->Bx)(ip->B_par,rr1,p);
  1160.       rr1 = p;
  1161.    }
  1162.    (ip->ATx)(ip->AT_par,rr1,r);        /* r = A^T*B*(b-A*x)  */
  1163.  
  1164.  
  1165.    old_inner = 0.0;
  1166.    for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
  1167.    {
  1168.       rr1 = r;
  1169.       if ( ip->Bx ) {
  1170.      (ip->Bx)(ip->B_par,r,z);        /* rr = B*r */
  1171.      rr1 = z;
  1172.       }
  1173.  
  1174.       inner = in_prod(r,rr1);
  1175.       nres = sqrt(fabs(inner));
  1176.       if (ip->info) ip->info(ip,nres,r,rr1);
  1177.       if ( ip->stop_crit(ip,nres,r,rr1) ) break;
  1178.  
  1179.       if ( ip->steps )    /* if ( ip->steps > 0 ) ... */
  1180.       {
  1181.      beta = inner/old_inner;
  1182.      p = v_mltadd(rr1,p,beta,p);
  1183.       }
  1184.       else        /* if ( ip->steps == 0 ) ... */
  1185.       {
  1186.      beta = 0.0;
  1187.      p = v_copy(rr1,p);
  1188.      old_inner = 0.0;
  1189.       }
  1190.       (ip->Ax)(ip->A_par,p,q);     /* q = A*p */
  1191.       if (ip->Bx) {
  1192.      (ip->Bx)(ip->B_par,q,z);
  1193.      (ip->ATx)(ip->AT_par,z,q);
  1194.      rr1 = q;            /* q = A^T*B*A*p */
  1195.       }
  1196.       else {
  1197.      (ip->ATx)(ip->AT_par,q,z);    /* z = A^T*A*p */
  1198.      rr1 = z;
  1199.       }
  1200.  
  1201.       alpha = inner/in_prod(rr1,p);
  1202.       v_mltadd(ip->x,p,alpha,ip->x);
  1203.       v_mltadd(r,rr1,-alpha,r);
  1204.       old_inner = inner;
  1205.    }
  1206.  
  1207.    return ip->x;
  1208. }
  1209.  
  1210. /* iter_spcgne -- a simple interface to iter_cgne() which 
  1211.    uses sparse matrix data structures
  1212.    -- assumes that B contains an actual preconditioner (or NULL)
  1213.    use always as follows:
  1214.       x = iter_spcgne(A,B,b,eps,x,limit,steps);
  1215.    or 
  1216.       x = iter_spcgne(A,B,b,eps,VNULL,limit,steps);
  1217.    In the second case the solution vector is created.
  1218. */
  1219. VEC  *iter_spcgne(A,B,b,eps,x,limit,steps)
  1220. SPMAT    *A, *B;
  1221. VEC    *b, *x;
  1222. double    eps;
  1223. int *steps, limit;
  1224. {    
  1225.    ITER *ip;
  1226.    
  1227.    ip = iter_get(0,0);
  1228.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  1229.    ip->A_par = (void *)A;
  1230.    ip->ATx = (Fun_Ax) sp_vm_mlt;
  1231.    ip->AT_par = (void *)A;
  1232.    if (B) {
  1233.       ip->Bx = (Fun_Ax) sp_mv_mlt;
  1234.       ip->B_par = (void *)B;
  1235.    }
  1236.    else {
  1237.       ip->Bx = (Fun_Ax) NULL;
  1238.       ip->B_par = NULL;
  1239.    }
  1240.    ip->info = (Fun_info) NULL;
  1241.    ip->b = b;
  1242.    ip->eps = eps;
  1243.    ip->limit = limit;
  1244.    ip->x = x;
  1245.    iter_cgne(ip);
  1246.    x = ip->x;
  1247.    if (steps) *steps = ip->steps;
  1248.    ip->shared_x = ip->shared_b = TRUE;
  1249.    iter_free(ip);   /* release only ITER structure */
  1250.    return x;        
  1251. }
  1252.  
  1253.  
  1254.  
  1255.  
  1256.