home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / itersym < prev    next >
Encoding:
Text File  |  1994-01-13  |  13.5 KB  |  585 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. /* itersym.c 17/09/93 */
  28.  
  29.  
  30. /* 
  31.   ITERATIVE METHODS - implementation of several iterative methods;
  32.   see also iter0.c
  33.   */
  34.  
  35. #include        <stdio.h>
  36. #include    <math.h>
  37. #include        "matrix.h"
  38. #include        "matrix2.h"
  39. #include    "sparse.h"
  40. #include        "iter.h"
  41.  
  42. static char rcsid[] = "$Id: itersym.c,v 1.1 1994/01/13 05:38:59 des Exp $";
  43.  
  44.  
  45. #ifdef ANSI_C
  46. VEC    *spCHsolve(SPMAT *,VEC *,VEC *);
  47. VEC    *trieig(VEC *,VEC *,MAT *);
  48. #else
  49. VEC    *spCHsolve();
  50. VEC    *trieig();
  51. #endif
  52.  
  53.  
  54.  
  55. /* iter_spcg -- a simple interface to iter_cg() which uses sparse matrix
  56.    data structures
  57.    -- assumes that LLT contains the Cholesky factorisation of the
  58.    actual preconditioner;
  59.    use always as follows:
  60.    x = iter_spcg(A,LLT,b,eps,x,limit,steps);
  61.    or 
  62.    x = iter_spcg(A,LLT,b,eps,VNULL,limit,steps);
  63.    In the second case the solution vector is created.
  64.    */
  65. VEC  *iter_spcg(A,LLT,b,eps,x,limit,steps)
  66. SPMAT    *A, *LLT;
  67. VEC    *b, *x;
  68. double    eps;
  69. int *steps, limit;
  70. {    
  71.    ITER *ip;
  72.    
  73.    ip = iter_get(0,0);
  74.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  75.    ip->A_par = (void *)A;
  76.    ip->Bx = (Fun_Ax) spCHsolve;
  77.    ip->B_par = (void *)LLT;
  78.    ip->info = (Fun_info) NULL;
  79.    ip->b = b;
  80.    ip->eps = eps;
  81.    ip->limit = limit;
  82.    ip->x = x;
  83.    iter_cg(ip);
  84.    x = ip->x;
  85.    if (steps) *steps = ip->steps;
  86.    ip->shared_x = ip->shared_b = TRUE;
  87.    iter_free(ip);   /* release only ITER structure */
  88.    return x;        
  89. }
  90.  
  91. /* 
  92.   Conjugate gradients method;
  93.   */
  94. VEC  *iter_cg(ip)
  95. ITER *ip;
  96. {
  97.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  98.    Real    alpha, beta, inner, old_inner, nres;
  99.    VEC *rr;   /* rr == r or rr == z */
  100.    
  101.    if (ip == INULL)
  102.      error(E_NULL,"iter_cg");
  103.    if (!ip->Ax || !ip->b)
  104.      error(E_NULL,"iter_cg");
  105.    if ( ip->x == ip->b )
  106.      error(E_INSITU,"iter_cg");
  107.    if (!ip->stop_crit)
  108.      error(E_NULL,"iter_cg");
  109.    
  110.    if ( ip->eps <= 0.0 )
  111.      ip->eps = MACHEPS;
  112.    
  113.    r = v_resize(r,ip->b->dim);
  114.    p = v_resize(p,ip->b->dim);
  115.    q = v_resize(q,ip->b->dim);
  116.    
  117.    MEM_STAT_REG(r,TYPE_VEC);
  118.    MEM_STAT_REG(p,TYPE_VEC);
  119.    MEM_STAT_REG(q,TYPE_VEC);
  120.    
  121.    if (ip->Bx != (Fun_Ax)NULL) {
  122.       z = v_resize(z,ip->b->dim);
  123.       MEM_STAT_REG(z,TYPE_VEC);
  124.       rr = z;
  125.    }
  126.    else rr = r;
  127.    
  128.    if (ip->x != VNULL) {
  129.       if (ip->x->dim != ip->b->dim)
  130.     error(E_SIZES,"iter_cg");
  131.       ip->Ax(ip->A_par,ip->x,p);            /* p = A*x */
  132.       v_sub(ip->b,p,r);                 /* r = b - A*x */
  133.    }
  134.    else {  /* ip->x == 0 */
  135.       ip->x = v_get(ip->b->dim);
  136.       ip->shared_x = FALSE;
  137.       v_copy(ip->b,r);
  138.    }
  139.    
  140.    old_inner = 0.0;
  141.    for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
  142.    {
  143.       if ( ip->Bx )
  144.     (ip->Bx)(ip->B_par,r,rr);        /* rr = B*r */
  145.       
  146.       inner = in_prod(rr,r);
  147.       nres = sqrt(fabs(inner));
  148.       if (ip->info) ip->info(ip,nres,r,rr);
  149.       if ( ip->stop_crit(ip,nres,r,rr) ) break;
  150.       
  151.       if ( ip->steps )    /* if ( ip->steps > 0 ) ... */
  152.       {
  153.      beta = inner/old_inner;
  154.      p = v_mltadd(rr,p,beta,p);
  155.       }
  156.       else        /* if ( ip->steps == 0 ) ... */
  157.       {
  158.      beta = 0.0;
  159.      p = v_copy(rr,p);
  160.      old_inner = 0.0;
  161.       }
  162.       (ip->Ax)(ip->A_par,p,q);     /* q = A*p */
  163.       alpha = inner/in_prod(p,q);
  164.       v_mltadd(ip->x,p,alpha,ip->x);
  165.       v_mltadd(r,q,-alpha,r);
  166.       old_inner = inner;
  167.    }
  168.    
  169.    return ip->x;
  170. }
  171.  
  172.  
  173.  
  174. /* iter_lanczos -- raw lanczos algorithm -- no re-orthogonalisation
  175.    -- creates T matrix of size == m,
  176.    but no larger than before beta_k == 0
  177.    -- uses passed routine to do matrix-vector multiplies */
  178. void    iter_lanczos(ip,a,b,beta2,Q)
  179. ITER    *ip;
  180. VEC    *a, *b;
  181. Real    *beta2;
  182. MAT    *Q;
  183. {
  184.    int    j;
  185.    static VEC    *v = VNULL, *w = VNULL, *tmp = VNULL;
  186.    Real    alpha, beta, c;
  187.    
  188.    if ( ! ip )
  189.      error(E_NULL,"iter_lanczos");
  190.    if ( ! ip->Ax || ! ip->x || ! a || ! b )
  191.      error(E_NULL,"iter_lanczos");
  192.    if ( ip->k <= 0 )
  193.      error(E_BOUNDS,"iter_lanczos");
  194.    if ( Q && ( Q->n < ip->x->dim || Q->m < ip->k ) )
  195.      error(E_SIZES,"iter_lanczos");
  196.    
  197.    a = v_resize(a,(u_int)ip->k);    
  198.    b = v_resize(b,(u_int)(ip->k-1));
  199.    v = v_resize(v,ip->x->dim);
  200.    w = v_resize(w,ip->x->dim);
  201.    tmp = v_resize(tmp,ip->x->dim);
  202.    MEM_STAT_REG(v,TYPE_VEC);
  203.    MEM_STAT_REG(w,TYPE_VEC);
  204.    MEM_STAT_REG(tmp,TYPE_VEC);
  205.    
  206.    beta = 1.0;
  207.    v_zero(a);
  208.    v_zero(b);
  209.    if (Q) m_zero(Q);
  210.    
  211.    /* normalise x as w */
  212.    c = v_norm2(ip->x);
  213.    if (c <= MACHEPS) { /* ip->x == 0 */
  214.       *beta2 = 0.0;
  215.       return;
  216.    }
  217.    else 
  218.      sv_mlt(1.0/c,ip->x,w);
  219.    
  220.    (ip->Ax)(ip->A_par,w,v);
  221.    
  222.    for ( j = 0; j < ip->k; j++ )
  223.    {
  224.       /* store w in Q if Q not NULL */
  225.       if ( Q ) set_row(Q,j,w);
  226.       
  227.       alpha = in_prod(w,v);
  228.       a->ve[j] = alpha;
  229.       v_mltadd(v,w,-alpha,v);
  230.       beta = v_norm2(v);
  231.       if ( beta == 0.0 )
  232.       {
  233.      *beta2 = 0.0;
  234.      return;
  235.       }
  236.       
  237.       if ( j < ip->k-1 )
  238.     b->ve[j] = beta;
  239.       v_copy(w,tmp);
  240.       sv_mlt(1/beta,v,w);
  241.       sv_mlt(-beta,tmp,v);
  242.       (ip->Ax)(ip->A_par,w,tmp);
  243.       v_add(v,tmp,v);
  244.    }
  245.    *beta2 = beta;
  246.    
  247. }
  248.  
  249. /* iter_splanczos -- version that uses sparse matrix data structure */
  250. void    iter_splanczos(A,m,x0,a,b,beta2,Q)
  251. SPMAT    *A;
  252. int     m;
  253. VEC     *x0, *a, *b;
  254. Real    *beta2;
  255. MAT     *Q;
  256. {    
  257.    ITER *ip;
  258.    
  259.    ip = iter_get(0,0);
  260.    ip->shared_x = ip->shared_b = TRUE;
  261.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  262.    ip->A_par = (void *) A;
  263.    ip->x = x0;
  264.    ip->k = m;
  265.    iter_lanczos(ip,a,b,beta2,Q);    
  266.    iter_free(ip);   /* release only ITER structure */
  267. }
  268.  
  269.  
  270.  
  271. extern    double    frexp(), ldexp();
  272.  
  273. /* product -- returns the product of a long list of numbers
  274.    -- answer stored in mant (mantissa) and expt (exponent) */
  275. static    double    product(a,offset,expt)
  276. VEC    *a;
  277. double    offset;
  278. int    *expt;
  279. {
  280.    Real    mant, tmp_fctr;
  281.    int    i, tmp_expt;
  282.    
  283.    if ( ! a )
  284.      error(E_NULL,"product");
  285.    
  286.    mant = 1.0;
  287.    *expt = 0;
  288.    if ( offset == 0.0 )
  289.      for ( i = 0; i < a->dim; i++ )
  290.      {
  291.     mant *= frexp(a->ve[i],&tmp_expt);
  292.     *expt += tmp_expt;
  293.     if ( ! (i % 10) )
  294.     {
  295.        mant = frexp(mant,&tmp_expt);
  296.        *expt += tmp_expt;
  297.     }
  298.      }
  299.    else
  300.      for ( i = 0; i < a->dim; i++ )
  301.      {
  302.     tmp_fctr = a->ve[i] - offset;
  303.     tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset :
  304.       MACHEPS*offset;
  305.     mant *= frexp(tmp_fctr,&tmp_expt);
  306.     *expt += tmp_expt;
  307.     if ( ! (i % 10) )
  308.     {
  309.        mant = frexp(mant,&tmp_expt);
  310.        *expt += tmp_expt;
  311.     }
  312.      }
  313.    
  314.    mant = frexp(mant,&tmp_expt);
  315.    *expt += tmp_expt;
  316.    
  317.    return mant;
  318. }
  319.  
  320. /* product2 -- returns the product of a long list of numbers
  321.    -- answer stored in mant (mantissa) and expt (exponent) */
  322. static    double    product2(a,k,expt)
  323. VEC    *a;
  324. int    k;    /* entry of a to leave out */
  325. int    *expt;
  326. {
  327.    Real    mant, mu, tmp_fctr;
  328.    int    i, tmp_expt;
  329.    
  330.    if ( ! a )
  331.      error(E_NULL,"product2");
  332.    if ( k < 0 || k >= a->dim )
  333.      error(E_BOUNDS,"product2");
  334.    
  335.    mant = 1.0;
  336.    *expt = 0;
  337.    mu = a->ve[k];
  338.    for ( i = 0; i < a->dim; i++ )
  339.    {
  340.       if ( i == k )
  341.     continue;
  342.       tmp_fctr = a->ve[i] - mu;
  343.       tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu;
  344.       mant *= frexp(tmp_fctr,&tmp_expt);
  345.       *expt += tmp_expt;
  346.       if ( ! (i % 10) )
  347.       {
  348.      mant = frexp(mant,&tmp_expt);
  349.      *expt += tmp_expt;
  350.       }
  351.    }
  352.    mant = frexp(mant,&tmp_expt);
  353.    *expt += tmp_expt;
  354.    
  355.    return mant;
  356. }
  357.  
  358. /* dbl_cmp -- comparison function to pass to qsort() */
  359. static    int    dbl_cmp(x,y)
  360. Real    *x, *y;
  361. {
  362.    Real    tmp;
  363.    
  364.    tmp = *x - *y;
  365.    return (tmp > 0 ? 1 : tmp < 0 ? -1: 0);
  366. }
  367.  
  368. /* iter_lanczos2 -- lanczos + error estimate for every e-val
  369.    -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978
  370.    -- returns multiple e-vals where multiple e-vals may not exist
  371.    -- returns evals vector */
  372. VEC    *iter_lanczos2(ip,evals,err_est)
  373. ITER     *ip;            /* ITER structure */
  374. VEC    *evals;        /* eigenvalue vector */
  375. VEC    *err_est;    /* error estimates of eigenvalues */
  376. {
  377.    VEC        *a;
  378.    static    VEC    *b=VNULL, *a2=VNULL, *b2=VNULL;
  379.    Real    beta, pb_mant, det_mant, det_mant1, det_mant2;
  380.    int    i, pb_expt, det_expt, det_expt1, det_expt2;
  381.    
  382.    if ( ! ip )
  383.      error(E_NULL,"iter_lanczos2");
  384.    if ( ! ip->Ax || ! ip->x )
  385.      error(E_NULL,"iter_lanczos2");
  386.    if ( ip->k <= 0 )
  387.      error(E_RANGE,"iter_lanczos2");
  388.    
  389.    a = evals;
  390.    a = v_resize(a,(u_int)ip->k);
  391.    b = v_resize(b,(u_int)(ip->k-1));
  392.    MEM_STAT_REG(b,TYPE_VEC);
  393.    
  394.    iter_lanczos(ip,a,b,&beta,MNULL);
  395.    
  396.    /* printf("# beta =%g\n",beta); */
  397.    pb_mant = 0.0;
  398.    if ( err_est )
  399.    {
  400.       pb_mant = product(b,(double)0.0,&pb_expt);
  401.       /* printf("# pb_mant = %g, pb_expt = %d\n",pb_mant, pb_expt); */
  402.    }
  403.    
  404.    /* printf("# diags =\n");    v_output(a); */
  405.    /* printf("# off diags =\n");    v_output(b); */
  406.    a2 = v_resize(a2,a->dim - 1);
  407.    b2 = v_resize(b2,b->dim - 1);
  408.    MEM_STAT_REG(a2,TYPE_VEC);
  409.    MEM_STAT_REG(b2,TYPE_VEC);
  410.    for ( i = 0; i < a2->dim - 1; i++ )
  411.    {
  412.       a2->ve[i] = a->ve[i+1];
  413.       b2->ve[i] = b->ve[i+1];
  414.    }
  415.    a2->ve[a2->dim-1] = a->ve[a2->dim];
  416.    
  417.    trieig(a,b,MNULL);
  418.    
  419.    /* sort evals as a courtesy */
  420.    qsort((void *)(a->ve),(int)(a->dim),sizeof(Real),(int (*)())dbl_cmp);
  421.    
  422.    /* error estimates */
  423.    if ( err_est )
  424.    {
  425.       err_est = v_resize(err_est,(u_int)ip->k);
  426.       
  427.       trieig(a2,b2,MNULL);
  428.       /* printf("# a =\n");    v_output(a); */
  429.       /* printf("# a2 =\n");    v_output(a2); */
  430.       
  431.       for ( i = 0; i < a->dim; i++ )
  432.       {
  433.      det_mant1 = product2(a,i,&det_expt1);
  434.      det_mant2 = product(a2,(double)a->ve[i],&det_expt2);
  435.      /* printf("# det_mant1=%g, det_expt1=%d\n",
  436.         det_mant1,det_expt1); */
  437.      /* printf("# det_mant2=%g, det_expt2=%d\n",
  438.         det_mant2,det_expt2); */
  439.      if ( det_mant1 == 0.0 )
  440.      {   /* multiple e-val of T */
  441.         err_est->ve[i] = 0.0;
  442.         continue;
  443.      }
  444.      else if ( det_mant2 == 0.0 )
  445.      {
  446.         err_est->ve[i] = HUGE;
  447.         continue;
  448.      }
  449.      if ( (det_expt1 + det_expt2) % 2 )
  450.        /* if odd... */
  451.        det_mant = sqrt(2.0*fabs(det_mant1*det_mant2));
  452.      else /* if even... */
  453.        det_mant = sqrt(fabs(det_mant1*det_mant2));
  454.      det_expt = (det_expt1+det_expt2)/2;
  455.      err_est->ve[i] = fabs(beta*
  456.                    ldexp(pb_mant/det_mant,pb_expt-det_expt));
  457.       }
  458.    }
  459.    
  460.    return a;
  461. }
  462.  
  463. /* iter_splanczos2 -- version of iter_lanczos2() that uses sparse matrix data
  464.    structure */
  465.  
  466. VEC    *iter_splanczos2(A,m,x0,evals,err_est)
  467. SPMAT    *A;
  468. int     m;
  469. VEC    *x0;        /* initial vector */
  470. VEC    *evals;        /* eigenvalue vector */
  471. VEC    *err_est;    /* error estimates of eigenvalues */
  472. {    
  473.    ITER *ip;
  474.    VEC *a;
  475.    
  476.    ip = iter_get(0,0);
  477.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  478.    ip->A_par = (void *) A;
  479.    ip->x = x0;
  480.    ip->k = m;
  481.    a = iter_lanczos2(ip,evals,err_est);    
  482.    ip->shared_x = ip->shared_b = TRUE;
  483.    iter_free(ip);   /* release only ITER structure */
  484.    return a;
  485. }
  486.  
  487.  
  488.  
  489.  
  490. /*
  491.   Conjugate gradient method
  492.   Another variant - mainly for testing
  493.   */
  494.  
  495. VEC  *iter_cg1(ip)
  496. ITER *ip;
  497. {
  498.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  499.    Real    alpha;
  500.    double inner,nres;
  501.    VEC *rr;   /* rr == r or rr == z */
  502.    
  503.    if (ip == INULL)
  504.      error(E_NULL,"iter_cg");
  505.    if (!ip->Ax || !ip->b)
  506.      error(E_NULL,"iter_cg");
  507.    if ( ip->x == ip->b )
  508.      error(E_INSITU,"iter_cg");
  509.    if (!ip->stop_crit)
  510.      error(E_NULL,"iter_cg");
  511.    
  512.    if ( ip->eps <= 0.0 )
  513.      ip->eps = MACHEPS;
  514.    
  515.    r = v_resize(r,ip->b->dim);
  516.    p = v_resize(p,ip->b->dim);
  517.    q = v_resize(q,ip->b->dim);
  518.    
  519.    MEM_STAT_REG(r,TYPE_VEC);
  520.    MEM_STAT_REG(p,TYPE_VEC);
  521.    MEM_STAT_REG(q,TYPE_VEC);
  522.    
  523.    if (ip->Bx != (Fun_Ax)NULL) {
  524.       z = v_resize(z,ip->b->dim);
  525.       MEM_STAT_REG(z,TYPE_VEC);
  526.       rr = z;
  527.    }
  528.    else rr = r;
  529.    
  530.    if (ip->x != VNULL) {
  531.       if (ip->x->dim != ip->b->dim)
  532.     error(E_SIZES,"iter_cg");
  533.       ip->Ax(ip->A_par,ip->x,p);            /* p = A*x */
  534.       v_sub(ip->b,p,r);                 /* r = b - A*x */
  535.    }
  536.    else {  /* ip->x == 0 */
  537.       ip->x = v_get(ip->b->dim);
  538.       ip->shared_x = FALSE;
  539.       v_copy(ip->b,r);
  540.    }
  541.    
  542.    if (ip->Bx) (ip->Bx)(ip->B_par,r,p);
  543.    else v_copy(r,p);
  544.    
  545.    inner = in_prod(p,r);
  546.    nres = sqrt(fabs(inner));
  547.    if (ip->info) ip->info(ip,nres,r,p);
  548.    if ( ip->stop_crit(ip,nres,r,p) ) return ip->x;
  549.    
  550.    for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
  551.    {
  552.       ip->Ax(ip->A_par,p,q);
  553.       inner = in_prod(q,p);
  554.       if (inner <= 0.0) {
  555.      warning(WARN_RES_LESS_0,"iter_cg");
  556.      break;
  557.       }
  558.       alpha = in_prod(p,r)/inner;
  559.       v_mltadd(ip->x,p,alpha,ip->x);
  560.       v_mltadd(r,q,-alpha,r);
  561.       
  562.       rr = r;
  563.       if (ip->Bx) {
  564.      ip->Bx(ip->B_par,r,z);
  565.      rr = z;
  566.       }
  567.       
  568.       nres = in_prod(r,rr);
  569.       if (nres < 0.0) {
  570.      warning(WARN_RES_LESS_0,"iter_cg");
  571.      break;
  572.       }
  573.       nres = sqrt(fabs(nres));
  574.       if (ip->info) ip->info(ip,nres,r,z);
  575.       if ( ip->stop_crit(ip,nres,r,z) ) break;
  576.       
  577.       alpha = -in_prod(rr,q)/inner;
  578.       v_mltadd(rr,p,alpha,p);
  579.       
  580.    }
  581.    
  582.    return ip->x;
  583. }
  584.  
  585.