home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / iternsym < prev    next >
Encoding:
Text File  |  1994-02-02  |  29.4 KB  |  1,237 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 Also the member ip->init_res is updated indirectly by 
  428.      ip->stop_crit.
  429.      */
  430.       if (ip->steps == 0) {                /* information for a user */
  431.      if (ip->info) (*ip->info)(ip,nres,As,rr); 
  432.      if ( (*ip->stop_crit)(ip,nres,As,rr) ) { 
  433.           /* iterative process is finished */
  434.         done = TRUE; 
  435.         break;
  436.      }
  437.       }
  438.       else if (nres <= 0.0) 
  439.     break;  /* residual is zero -> finish the process */ 
  440.       
  441.       /* save this residual in the first row of N */
  442.       v.ve = N->me[0];
  443.       v_copy(rr,&v);
  444.       
  445.       for (i = 0; i < ip->k && ip->steps <= ip->limit; i++) {
  446.      ip->steps++;
  447.      v.ve = N->me[i];                /* pointer to a row of N (=s_i) */
  448.      /* note that we must use here &v, not v */
  449.      (*ip->Ax)(ip->A_par,&v,As); 
  450.      rr = As;                        /* As = A*s_i */
  451.      if (ip->Bx) {
  452.         (*ip->Bx)(ip->B_par,As,z);    /* z = B*A*s_i  */
  453.         rr = z;
  454.      }
  455.      
  456.      if (i < ip->k - 1) {
  457.         s.ve = N->me[i+1];         /* pointer to a row of N (=s_{i+1}) */
  458.         v_copy(rr,&s);                   /* s_{i+1} = B*A*s_i */
  459.         for (j = 0; j <= i-1; j++) {
  460.            v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
  461.            beta->ve[j] = in_prod(&v,rr);    /* beta_{j,i} */
  462.                    /* s_{i+1} -= beta_{j,i}*s_{j+1} */
  463.            v_mltadd(&s,&v,- beta->ve[j],&s);    
  464.         }
  465.         
  466.          /* beta_{i,i} = ||s_{i+1}||_2 */
  467.         beta->ve[i] = nres = v_norm2(&s);     
  468.         if ( nres <= 0.0) break;         /* if s_{i+1} == 0 then finish */
  469.         sv_mlt(1.0/nres,&s,&s);           /* normalize s_{i+1} */
  470.         
  471.         v.ve = N->me[0];
  472.         alpha->ve[i] = in_prod(&v,&s);     /* alpha_i = (s_0 , s_{i+1}) */
  473.         
  474.      }
  475.      else {
  476.         for (j = 0; j <= i-1; j++) {
  477.            v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
  478.            beta->ve[j] = in_prod(&v,rr);       /* beta_{j,i} */
  479.         }
  480.         
  481.         nres = in_prod(rr,rr);                 /* rr = B*A*s_{k-1} */
  482.         for (j = 0; j <= i-1; j++)
  483.               nres -= beta->ve[j]*beta->ve[j];
  484.         if (nres <= 0.0)  break;               /* s_k is zero */
  485.         else beta->ve[i] = sqrt(nres);         /* beta_{k-1,k-1} */
  486.         
  487.         v.ve = N->me[0];
  488.         alpha->ve[i] = in_prod(&v,rr); 
  489.         for (j = 0; j <= i-1; j++)
  490.               alpha->ve[i] -= beta->ve[j]*alpha->ve[j];
  491.         alpha->ve[i] /= beta->ve[i];                /* alpha_{k-1} */
  492.         
  493.      }
  494.      
  495.      set_col(H,i,beta);
  496.      
  497.      dd -= alpha->ve[i]*alpha->ve[i];
  498.      nres = sqrt(fabs((double) dd));
  499.      if (dd < 0.0)  {     /* do restart */
  500.         if (ip->info) (*ip->info)(ip,-nres,VNULL,VNULL);  
  501.         break;
  502.      }
  503.      
  504.      if (ip->info) (*ip->info)(ip,nres,VNULL,VNULL);     
  505.      if ( ip->stop_crit(ip,nres,VNULL,VNULL) ) {
  506.         /* stopping criterion is satisfied */
  507.         done = TRUE;
  508.         break;
  509.      }
  510.      
  511.       } /* end of for */
  512.       
  513.       if (nres <= 0.0) {
  514.      i--;
  515.      done = TRUE;
  516.       }
  517.       if (i >= ip->k) i = ip->k - 1;
  518.       
  519.       /* use (i+1) by (i+1) submatrix of H */
  520.       H = m_resize(H,i+1,i+1);
  521.       alpha = v_resize(alpha,i+1);
  522.       Usolve(H,alpha,alpha,0.0);       /* c_i is saved in alpha */
  523.       
  524.       for (j = 0; j <= i; j++) {
  525.      v.ve = N->me[j];
  526.      v_mltadd(ip->x,&v,alpha->ve[j],ip->x);
  527.       }
  528.       
  529.       
  530.       if (done) break;              /* stop the iterative process */
  531.       alpha = v_resize(alpha,ip->k);
  532.       H = m_resize(H,ip->k,ip->k);
  533.       
  534.    }  /* end of while */
  535.    
  536.    return ip->x;                    /* return the solution */
  537. }
  538.  
  539.  
  540.  
  541. /* iter_spmgcr - a simple interface to iter_mgcr */
  542. /* no preconditioner */
  543. VEC    *iter_spmgcr(A,B,b,tol,x,k,limit,steps)
  544. SPMAT    *A, *B;
  545. VEC    *b, *x;
  546. double    tol;
  547. int *steps,k,limit;
  548. {
  549.    ITER *ip;
  550.    
  551.    ip = iter_get(0,0);
  552.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  553.    ip->A_par = (void *) A;
  554.    if (B) {
  555.       ip->Bx = (Fun_Ax) sp_mv_mlt;
  556.       ip->B_par = (void *) B;
  557.    }
  558.    else {
  559.       ip->Bx = (Fun_Ax) NULL;
  560.       ip->B_par = NULL;
  561.    }
  562.  
  563.    ip->k = k;
  564.    ip->limit = limit;
  565.    ip->info = (Fun_info) NULL;
  566.    ip->b = b;
  567.    ip->eps = tol;
  568.    ip->x = x;
  569.    iter_mgcr(ip);
  570.    x = ip->x;
  571.    if (steps) *steps = ip->steps;
  572.    ip->shared_x = ip->shared_b = TRUE;
  573.    iter_free(ip);   /* release only ITER structure */
  574.    return x;        
  575. }
  576.  
  577.  
  578.  
  579. /* 
  580.   Conjugate gradients method for a normal equation
  581.   a preconditioner B must be symmetric !!
  582. */
  583. VEC  *iter_cgne(ip)
  584. ITER *ip;
  585. {
  586.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  587.    Real    alpha, beta, inner, old_inner, nres;
  588.    VEC *rr1;   /* pointer only */
  589.    
  590.    if (ip == INULL)
  591.      error(E_NULL,"iter_cgne");
  592.    if (!ip->Ax || ! ip->ATx || !ip->b)
  593.      error(E_NULL,"iter_cgne");
  594.    if ( ip->x == ip->b )
  595.      error(E_INSITU,"iter_cgne");
  596.    if (!ip->stop_crit)
  597.      error(E_NULL,"iter_cgne");
  598.    
  599.    if ( ip->eps <= 0.0 )
  600.      ip->eps = MACHEPS;
  601.    
  602.    r = v_resize(r,ip->b->dim);
  603.    p = v_resize(p,ip->b->dim);
  604.    q = v_resize(q,ip->b->dim);
  605.  
  606.    MEM_STAT_REG(r,TYPE_VEC);
  607.    MEM_STAT_REG(p,TYPE_VEC);
  608.    MEM_STAT_REG(q,TYPE_VEC);
  609.  
  610.    z = v_resize(z,ip->b->dim);
  611.    MEM_STAT_REG(z,TYPE_VEC);
  612.  
  613.    if (ip->x) {
  614.       if (ip->x->dim != ip->b->dim)
  615.     error(E_SIZES,"iter_cgne");
  616.       ip->Ax(ip->A_par,ip->x,p);            /* p = A*x */
  617.       v_sub(ip->b,p,z);                 /* z = b - A*x */
  618.    }
  619.    else {  /* ip->x == 0 */
  620.       ip->x = v_get(ip->b->dim);
  621.       ip->shared_x = FALSE;
  622.       v_copy(ip->b,z);
  623.    }
  624.    rr1 = z;
  625.    if (ip->Bx) {
  626.       (ip->Bx)(ip->B_par,rr1,p);
  627.       rr1 = p;
  628.    }
  629.    (ip->ATx)(ip->AT_par,rr1,r);        /* r = A^T*B*(b-A*x)  */
  630.  
  631.  
  632.    old_inner = 0.0;
  633.    for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
  634.    {
  635.       rr1 = r;
  636.       if ( ip->Bx ) {
  637.      (ip->Bx)(ip->B_par,r,z);        /* rr = B*r */
  638.      rr1 = z;
  639.       }
  640.  
  641.       inner = in_prod(r,rr1);
  642.       nres = sqrt(fabs(inner));
  643.       if (ip->info) ip->info(ip,nres,r,rr1);
  644.       if ( ip->stop_crit(ip,nres,r,rr1) ) break;
  645.  
  646.       if ( ip->steps )    /* if ( ip->steps > 0 ) ... */
  647.       {
  648.      beta = inner/old_inner;
  649.      p = v_mltadd(rr1,p,beta,p);
  650.       }
  651.       else        /* if ( ip->steps == 0 ) ... */
  652.       {
  653.      beta = 0.0;
  654.      p = v_copy(rr1,p);
  655.      old_inner = 0.0;
  656.       }
  657.       (ip->Ax)(ip->A_par,p,q);     /* q = A*p */
  658.       if (ip->Bx) {
  659.      (ip->Bx)(ip->B_par,q,z);
  660.      (ip->ATx)(ip->AT_par,z,q);
  661.      rr1 = q;            /* q = A^T*B*A*p */
  662.       }
  663.       else {
  664.      (ip->ATx)(ip->AT_par,q,z);    /* z = A^T*A*p */
  665.      rr1 = z;
  666.       }
  667.  
  668.       alpha = inner/in_prod(rr1,p);
  669.       v_mltadd(ip->x,p,alpha,ip->x);
  670.       v_mltadd(r,rr1,-alpha,r);
  671.       old_inner = inner;
  672.    }
  673.  
  674.    return ip->x;
  675. }
  676.  
  677. /* iter_spcgne -- a simple interface to iter_cgne() which 
  678.    uses sparse matrix data structures
  679.    -- assumes that B contains an actual preconditioner (or NULL)
  680.    use always as follows:
  681.       x = iter_spcgne(A,B,b,eps,x,limit,steps);
  682.    or 
  683.       x = iter_spcgne(A,B,b,eps,VNULL,limit,steps);
  684.    In the second case the solution vector is created.
  685. */
  686. VEC  *iter_spcgne(A,B,b,eps,x,limit,steps)
  687. SPMAT    *A, *B;
  688. VEC    *b, *x;
  689. double    eps;
  690. int *steps, limit;
  691. {    
  692.    ITER *ip;
  693.    
  694.    ip = iter_get(0,0);
  695.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  696.    ip->A_par = (void *)A;
  697.    ip->ATx = (Fun_Ax) sp_vm_mlt;
  698.    ip->AT_par = (void *)A;
  699.    if (B) {
  700.       ip->Bx = (Fun_Ax) sp_mv_mlt;
  701.       ip->B_par = (void *)B;
  702.    }
  703.    else {
  704.       ip->Bx = (Fun_Ax) NULL;
  705.       ip->B_par = NULL;
  706.    }
  707.    ip->info = (Fun_info) NULL;
  708.    ip->b = b;
  709.    ip->eps = eps;
  710.    ip->limit = limit;
  711.    ip->x = x;
  712.    iter_cgne(ip);
  713.    x = ip->x;
  714.    if (steps) *steps = ip->steps;
  715.    ip->shared_x = ip->shared_b = TRUE;
  716.    iter_free(ip);   /* release only ITER structure */
  717.    return x;        
  718. }
  719.  
  720.  
  721.  
  722.  
  723. FileDataŵitersym6Eÿÿÿ8S@V´
  724. /**************************************************************************
  725. **
  726. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  727. **
  728. **                 Meschach Library
  729. ** 
  730. ** This Meschach Library is provided "as is" without any express 
  731. ** or implied warranty of any kind with respect to this software. 
  732. ** In particular the authors shall not be liable for any direct, 
  733. ** indirect, special, incidental or consequential damages arising 
  734. ** in any way from use of the software.
  735. ** 
  736. ** Everyone is granted permission to copy, modify and redistribute this
  737. ** Meschach Library, provided:
  738. **  1.  All copies contain this copyright notice.
  739. **  2.  All modified copies shall carry a notice stating who
  740. **      made the last modification and the date of such modification.
  741. **  3.  No charge is made for this software or works derived from it.  
  742. **      This clause shall not be construed as constraining other software
  743. **      distributed on the same medium as this software, nor is a
  744. **      distribution fee considered a charge.
  745. **
  746. ***************************************************************************/
  747.  
  748.  
  749. /* itersym.c 17/09/93 */
  750.  
  751.  
  752. /* 
  753.   ITERATIVE METHODS - implementation of several iterative methods;
  754.   see also iter0.c
  755.   */
  756.  
  757. #include        <stdio.h>
  758. #include    <math.h>
  759. #include        "matrix.h"
  760. #include        "matrix2.h"
  761. #include    "sparse.h"
  762. #include        "iter.h"
  763.  
  764. static char rcsid[] = "$Id: itersym.c,v 1.1 1994/01/13 05:38:59 des Exp $";
  765.  
  766.  
  767. #ifdef ANSI_C
  768. VEC    *spCHsolve(SPMAT *,VEC *,VEC *);
  769. VEC    *trieig(VEC *,VEC *,MAT *);
  770. #else
  771. VEC    *spCHsolve();
  772. VEC    *trieig();
  773. #endif
  774.  
  775.  
  776.  
  777. /* iter_spcg -- a simple interface to iter_cg() which uses sparse matrix
  778.    data structures
  779.    -- assumes that LLT contains the Cholesky factorisation of the
  780.    actual preconditioner;
  781.    use always as follows:
  782.    x = iter_spcg(A,LLT,b,eps,x,limit,steps);
  783.    or 
  784.    x = iter_spcg(A,LLT,b,eps,VNULL,limit,steps);
  785.    In the second case the solution vector is created.
  786.    */
  787. VEC  *iter_spcg(A,LLT,b,eps,x,limit,steps)
  788. SPMAT    *A, *LLT;
  789. VEC    *b, *x;
  790. double    eps;
  791. int *steps, limit;
  792. {    
  793.    ITER *ip;
  794.    
  795.    ip = iter_get(0,0);
  796.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  797.    ip->A_par = (void *)A;
  798.    ip->Bx = (Fun_Ax) spCHsolve;
  799.    ip->B_par = (void *)LLT;
  800.    ip->info = (Fun_info) NULL;
  801.    ip->b = b;
  802.    ip->eps = eps;
  803.    ip->limit = limit;
  804.    ip->x = x;
  805.    iter_cg(ip);
  806.    x = ip->x;
  807.    if (steps) *steps = ip->steps;
  808.    ip->shared_x = ip->shared_b = TRUE;
  809.    iter_free(ip);   /* release only ITER structure */
  810.    return x;        
  811. }
  812.  
  813. /* 
  814.   Conjugate gradients method;
  815.   */
  816. VEC  *iter_cg(ip)
  817. ITER *ip;
  818. {
  819.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  820.    Real    alpha, beta, inner, old_inner, nres;
  821.    VEC *rr;   /* rr == r or rr == z */
  822.    
  823.    if (ip == INULL)
  824.      error(E_NULL,"iter_cg");
  825.    if (!ip->Ax || !ip->b)
  826.      error(E_NULL,"iter_cg");
  827.    if ( ip->x == ip->b )
  828.      error(E_INSITU,"iter_cg");
  829.    if (!ip->stop_crit)
  830.      error(E_NULL,"iter_cg");
  831.    
  832.    if ( ip->eps <= 0.0 )
  833.      ip->eps = MACHEPS;
  834.    
  835.    r = v_resize(r,ip->b->dim);
  836.    p = v_resize(p,ip->b->dim);
  837.    q = v_resize(q,ip->b->dim);
  838.    
  839.    MEM_STAT_REG(r,TYPE_VEC);
  840.    MEM_STAT_REG(p,TYPE_VEC);
  841.    MEM_STAT_REG(q,TYPE_VEC);
  842.    
  843.    if (ip->Bx != (Fun_Ax)NULL) {
  844.       z = v_resize(z,ip->b->dim);
  845.       MEM_STAT_REG(z,TYPE_VEC);
  846.       rr = z;
  847.    }
  848.    else rr = r;
  849.    
  850.    if (ip->x != VNULL) {
  851.       if (ip->x->dim != ip->b->dim)
  852.     error(E_SIZES,"iter_cg");
  853.       ip->Ax(ip->A_par,ip->x,p);            /* p = A*x */
  854.       v_sub(ip->b,p,r);                 /* r = b - A*x */
  855.    }
  856.    else {  /* ip->x == 0 */
  857.       ip->x = v_get(ip->b->dim);
  858.       ip->shared_x = FALSE;
  859.       v_copy(ip->b,r);
  860.    }
  861.    
  862.    old_inner = 0.0;
  863.    for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
  864.    {
  865.       if ( ip->Bx )
  866.     (ip->Bx)(ip->B_par,r,rr);        /* rr = B*r */
  867.       
  868.       inner = in_prod(rr,r);
  869.       nres = sqrt(fabs(inner));
  870.       if (ip->info) ip->info(ip,nres,r,rr);
  871.       if ( ip->stop_crit(ip,nres,r,rr) ) break;
  872.       
  873.       if ( ip->steps )    /* if ( ip->steps > 0 ) ... */
  874.       {
  875.      beta = inner/old_inner;
  876.      p = v_mltadd(rr,p,beta,p);
  877.       }
  878.       else        /* if ( ip->steps == 0 ) ... */
  879.       {
  880.      beta = 0.0;
  881.      p = v_copy(rr,p);
  882.      old_inner = 0.0;
  883.       }
  884.       (ip->Ax)(ip->A_par,p,q);     /* q = A*p */
  885.       alpha = inner/in_prod(p,q);
  886.       v_mltadd(ip->x,p,alpha,ip->x);
  887.       v_mltadd(r,q,-alpha,r);
  888.       old_inner = inner;
  889.    }
  890.    
  891.    return ip->x;
  892. }
  893.  
  894.  
  895.  
  896. /* iter_lanczos -- raw lanczos algorithm -- no re-orthogonalisation
  897.    -- creates T matrix of size == m,
  898.    but no larger than before beta_k == 0
  899.    -- uses passed routine to do matrix-vector multiplies */
  900. void    iter_lanczos(ip,a,b,beta2,Q)
  901. ITER    *ip;
  902. VEC    *a, *b;
  903. Real    *beta2;
  904. MAT    *Q;
  905. {
  906.    int    j;
  907.    static VEC    *v = VNULL, *w = VNULL, *tmp = VNULL;
  908.    Real    alpha, beta, c;
  909.    
  910.    if ( ! ip )
  911.      error(E_NULL,"iter_lanczos");
  912.    if ( ! ip->Ax || ! ip->x || ! a || ! b )
  913.      error(E_NULL,"iter_lanczos");
  914.    if ( ip->k <= 0 )
  915.      error(E_BOUNDS,"iter_lanczos");
  916.    if ( Q && ( Q->n < ip->x->dim || Q->m < ip->k ) )
  917.      error(E_SIZES,"iter_lanczos");
  918.    
  919.    a = v_resize(a,(u_int)ip->k);    
  920.    b = v_resize(b,(u_int)(ip->k-1));
  921.    v = v_resize(v,ip->x->dim);
  922.    w = v_resize(w,ip->x->dim);
  923.    tmp = v_resize(tmp,ip->x->dim);
  924.    MEM_STAT_REG(v,TYPE_VEC);
  925.    MEM_STAT_REG(w,TYPE_VEC);
  926.    MEM_STAT_REG(tmp,TYPE_VEC);
  927.    
  928.    beta = 1.0;
  929.    v_zero(a);
  930.    v_zero(b);
  931.    if (Q) m_zero(Q);
  932.    
  933.    /* normalise x as w */
  934.    c = v_norm2(ip->x);
  935.    if (c <= MACHEPS) { /* ip->x == 0 */
  936.       *beta2 = 0.0;
  937.       return;
  938.    }
  939.    else 
  940.      sv_mlt(1.0/c,ip->x,w);
  941.    
  942.    (ip->Ax)(ip->A_par,w,v);
  943.    
  944.    for ( j = 0; j < ip->k; j++ )
  945.    {
  946.       /* store w in Q if Q not NULL */
  947.       if ( Q ) set_row(Q,j,w);
  948.       
  949.       alpha = in_prod(w,v);
  950.       a->ve[j] = alpha;
  951.       v_mltadd(v,w,-alpha,v);
  952.       beta = v_norm2(v);
  953.       if ( beta == 0.0 )
  954.       {
  955.      *beta2 = 0.0;
  956.      return;
  957.       }
  958.       
  959.       if ( j < ip->k-1 )
  960.     b->ve[j] = beta;
  961.       v_copy(w,tmp);
  962.       sv_mlt(1/beta,v,w);
  963.       sv_mlt(-beta,tmp,v);
  964.       (ip->Ax)(ip->A_par,w,tmp);
  965.       v_add(v,tmp,v);
  966.    }
  967.    *beta2 = beta;
  968.    
  969. }
  970.  
  971. /* iter_splanczos -- version that uses sparse matrix data structure */
  972. void    iter_splanczos(A,m,x0,a,b,beta2,Q)
  973. SPMAT    *A;
  974. int     m;
  975. VEC     *x0, *a, *b;
  976. Real    *beta2;
  977. MAT     *Q;
  978. {    
  979.    ITER *ip;
  980.    
  981.    ip = iter_get(0,0);
  982.    ip->shared_x = ip->shared_b = TRUE;
  983.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  984.    ip->A_par = (void *) A;
  985.    ip->x = x0;
  986.    ip->k = m;
  987.    iter_lanczos(ip,a,b,beta2,Q);    
  988.    iter_free(ip);   /* release only ITER structure */
  989. }
  990.  
  991.  
  992.  
  993. extern    double    frexp(), ldexp();
  994.  
  995. /* product -- returns the product of a long list of numbers
  996.    -- answer stored in mant (mantissa) and expt (exponent) */
  997. static    double    product(a,offset,expt)
  998. VEC    *a;
  999. double    offset;
  1000. int    *expt;
  1001. {
  1002.    Real    mant, tmp_fctr;
  1003.    int    i, tmp_expt;
  1004.    
  1005.    if ( ! a )
  1006.      error(E_NULL,"product");
  1007.    
  1008.    mant = 1.0;
  1009.    *expt = 0;
  1010.    if ( offset == 0.0 )
  1011.      for ( i = 0; i < a->dim; i++ )
  1012.      {
  1013.     mant *= frexp(a->ve[i],&tmp_expt);
  1014.     *expt += tmp_expt;
  1015.     if ( ! (i % 10) )
  1016.     {
  1017.        mant = frexp(mant,&tmp_expt);
  1018.        *expt += tmp_expt;
  1019.     }
  1020.      }
  1021.    else
  1022.      for ( i = 0; i < a->dim; i++ )
  1023.      {
  1024.     tmp_fctr = a->ve[i] - offset;
  1025.     tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset :
  1026.       MACHEPS*offset;
  1027.     mant *= frexp(tmp_fctr,&tmp_expt);
  1028.     *expt += tmp_expt;
  1029.     if ( ! (i % 10) )
  1030.     {
  1031.        mant = frexp(mant,&tmp_expt);
  1032.        *expt += tmp_expt;
  1033.     }
  1034.      }
  1035.    
  1036.    mant = frexp(mant,&tmp_expt);
  1037.    *expt += tmp_expt;
  1038.    
  1039.    return mant;
  1040. }
  1041.  
  1042. /* product2 -- returns the product of a long list of numbers
  1043.    -- answer stored in mant (mantissa) and expt (exponent) */
  1044. static    double    product2(a,k,expt)
  1045. VEC    *a;
  1046. int    k;    /* entry of a to leave out */
  1047. int    *expt;
  1048. {
  1049.    Real    mant, mu, tmp_fctr;
  1050.    int    i, tmp_expt;
  1051.    
  1052.    if ( ! a )
  1053.      error(E_NULL,"product2");
  1054.    if ( k < 0 || k >= a->dim )
  1055.      error(E_BOUNDS,"product2");
  1056.    
  1057.    mant = 1.0;
  1058.    *expt = 0;
  1059.    mu = a->ve[k];
  1060.    for ( i = 0; i < a->dim; i++ )
  1061.    {
  1062.       if ( i == k )
  1063.     continue;
  1064.       tmp_fctr = a->ve[i] - mu;
  1065.       tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu;
  1066.       mant *= frexp(tmp_fctr,&tmp_expt);
  1067.       *expt += tmp_expt;
  1068.       if ( ! (i % 10) )
  1069.       {
  1070.      mant = frexp(mant,&tmp_expt);
  1071.      *expt += tmp_expt;
  1072.       }
  1073.    }
  1074.    mant = frexp(mant,&tmp_expt);
  1075.    *expt += tmp_expt;
  1076.    
  1077.    return mant;
  1078. }
  1079.  
  1080. /* dbl_cmp -- comparison function to pass to qsort() */
  1081. static    int    dbl_cmp(x,y)
  1082. Real    *x, *y;
  1083. {
  1084.    Real    tmp;
  1085.    
  1086.    tmp = *x - *y;
  1087.    return (tmp > 0 ? 1 : tmp < 0 ? -1: 0);
  1088. }
  1089.  
  1090. /* iter_lanczos2 -- lanczos + error estimate for every e-val
  1091.    -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978
  1092.    -- returns multiple e-vals where multiple e-vals may not exist
  1093.    -- returns evals vector */
  1094. VEC    *iter_lanczos2(ip,evals,err_est)
  1095. ITER     *ip;            /* ITER structure */
  1096. VEC    *evals;        /* eigenvalue vector */
  1097. VEC    *err_est;    /* error estimates of eigenvalues */
  1098. {
  1099.    VEC        *a;
  1100.    static    VEC    *b=VNULL, *a2=VNULL, *b2=VNULL;
  1101.    Real    beta, pb_mant, det_mant, det_mant1, det_mant2;
  1102.    int    i, pb_expt, det_expt, det_expt1, det_expt2;
  1103.    
  1104.    if ( ! ip )
  1105.      error(E_NULL,"iter_lanczos2");
  1106.    if ( ! ip->Ax || ! ip->x )
  1107.      error(E_NULL,"iter_lanczos2");
  1108.    if ( ip->k <= 0 )
  1109.      error(E_RANGE,"iter_lanczos2");
  1110.    
  1111.    a = evals;
  1112.    a = v_resize(a,(u_int)ip->k);
  1113.    b = v_resize(b,(u_int)(ip->k-1));
  1114.    MEM_STAT_REG(b,TYPE_VEC);
  1115.    
  1116.    iter_lanczos(ip,a,b,&beta,MNULL);
  1117.    
  1118.    /* printf("# beta =%g\n",beta); */
  1119.    pb_mant = 0.0;
  1120.    if ( err_est )
  1121.    {
  1122.       pb_mant = product(b,(double)0.0,&pb_expt);
  1123.       /* printf("# pb_mant = %g, pb_expt = %d\n",pb_mant, pb_expt); */
  1124.    }
  1125.    
  1126.    /* printf("# diags =\n");    v_output(a); */
  1127.    /* printf("# off diags =\n");    v_output(b); */
  1128.    a2 = v_resize(a2,a->dim - 1);
  1129.    b2 = v_resize(b2,b->dim - 1);
  1130.    MEM_STAT_REG(a2,TYPE_VEC);
  1131.    MEM_STAT_REG(b2,TYPE_VEC);
  1132.    for ( i = 0; i < a2->dim - 1; i++ )
  1133.    {
  1134.       a2->ve[i] = a->ve[i+1];
  1135.       b2->ve[i] = b->ve[i+1];
  1136.    }
  1137.    a2->ve[a2->dim-1] = a->ve[a2->dim];
  1138.    
  1139.    trieig(a,b,MNULL);
  1140.    
  1141.    /* sort evals as a courtesy */
  1142.    qsort((void *)(a->ve),(int)(a->dim),sizeof(Real),(int (*)())dbl_cmp);
  1143.    
  1144.    /* error estimates */
  1145.    if ( err_est )
  1146.    {
  1147.       err_est = v_resize(err_est,(u_int)ip->k);
  1148.       
  1149.       trieig(a2,b2,MNULL);
  1150.       /* printf("# a =\n");    v_output(a); */
  1151.       /* printf("# a2 =\n");    v_output(a2); */
  1152.       
  1153.       for ( i = 0; i < a->dim; i++ )
  1154.       {
  1155.      det_mant1 = product2(a,i,&det_expt1);
  1156.      det_mant2 = product(a2,(double)a->ve[i],&det_expt2);
  1157.      /*0  ,k+1,tmp1,v_entry(beta,k));
  1158.         /* printf("A = ");        m_output(A); */
  1159.     }
  1160.  
  1161.     return (A);
  1162. }
  1163.  
  1164. /* makeHQ -- construct the Hessenberg orthogonalising matrix Q;
  1165.     -- i.e. Hess M = Q.M.Q'    */
  1166. MAT    *makeHQ(H, diag, beta, Qout)
  1167. MAT    *H, *Qout;
  1168. VEC    *diag, *beta;
  1169. {
  1170.     int    i, j, limit;
  1171.     static    VEC    *tmp1 = VNULL, *tmp2 = VNULL;
  1172.  
  1173.     if ( H==(MAT *)NULL || diag==(VEC *)NULL || beta==(VEC *)NULL )
  1174.         error(E_NULL,"makeHQ");
  1175.     limit = H->m - 1;
  1176.     if ( diag->dim < limit || beta->dim < limit )
  1177.         error(E_SIZES,"makeHQ");
  1178.     if ( H->m != H->n )
  1179.         error(E_SQUARE,"makeHQ");
  1180.     Qout = m_resize(Qout,H->m,H->m);
  1181.  
  1182.     tmp1 = v_resize(tmp1,H->m);
  1183.     tmp2 = v_resize(tmp2,H->m);
  1184.     MEM_STAT_REG(tmp1,TYPE_VEC);
  1185.     MEM_STAT_REG(tmp2,TYPE_VEC);
  1186.  
  1187.     for ( i = 0; i < H->m; i++ )
  1188.     {
  1189.         /* tmp1 = i'th basis vector */
  1190.         for ( j = 0; j < H->m; j++ )
  1191.             /* tmp1->ve[j] = 0.0; */
  1192.             v_set_val(tmp1,j,0.0);
  1193.         /* tmp1->ve[i] = 1.0; */
  1194.         v_set_val(tmp1,i,1.0);
  1195.  
  1196.         /* apply H/h transforms in reverse order */
  1197.         for ( j = limit-1; j >= 0; j-- )
  1198.         {
  1199.             get_col(H,(u_int)j,tmp2);
  1200.             /* tmp2->ve[j+1] = diag->ve[j]; */
  1201.             v_set_val(tmp2,j+1,v_entry(diag,j));
  1202.             hhtrvec(tmp2,beta->ve[j],j+1,tmp1,tmp1);
  1203.         }
  1204.  
  1205.         /* insert into Qout */
  1206.         set_col(Qout,(u_int)i,tmp1);
  1207.     }
  1208.  
  1209.     return (Qout);
  1210. }
  1211.  
  1212. /* makeH -- construct actual Hessenberg matrix */
  1213. MAT    *makeH(H,Hout)
  1214. MAT    *H, *Hout;
  1215. {
  1216.     int    i, j, limit;
  1217.  
  1218.     if ( H==(MAT *)NULL )
  1219.         error(E_NULL,"makeH");
  1220.     if ( H->m != H->n )
  1221.         error(E_SQUARE,"makeH");
  1222.     Hout = m_resize(Hout,H->m,H->m);
  1223.     Hout = m_copy(H,Hout);
  1224.  
  1225.     limit = H->m;
  1226.     for ( i = 1; i < limit; i++ )
  1227.         for ( j = 0; j < i-1; j++ )
  1228.             /* Hout->me[i][j] = 0.0;*/
  1229.             m_set_val(Hout,i,j,0.0);
  1230.  
  1231.     return (Hout);
  1232. }
  1233.  
  1234. FileDataŵhsehldrEÿÿÿä@ï²
  1235. /**************************************************************************
  1236. **
  1237. ** Copyright (C) 19