home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / conjgrad < prev    next >
Encoding:
Text File  |  1994-01-13  |  8.1 KB  |  350 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. /*
  28.     Conjugate gradient routines file
  29.     Uses sparse matrix input & sparse Cholesky factorisation in pccg().
  30.  
  31.     All the following routines use routines to define a matrix
  32.         rather than use any explicit representation
  33.         (with the exeception of the pccg() pre-conditioner)
  34.     The matrix A is defined by
  35.  
  36.         VEC *(*A)(void *params, VEC *x, VEC *y)
  37.  
  38.     where y = A.x on exit, and y is returned. The params argument is
  39.     intended to make it easier to re-use & modify such routines.
  40.  
  41.     If we have a sparse matrix data structure
  42.         SPMAT    *A_mat;
  43.     then these can be used by passing sp_mv_mlt as the function, and
  44.     A_mat as the param.
  45. */
  46.  
  47. #include    <stdio.h>
  48. #include    <math.h>
  49. #include    "matrix.h"
  50. #include    "sparse.h"
  51. static char    rcsid[] = "$Id: conjgrad.c,v 1.4 1994/01/13 05:36:45 des Exp $";
  52.  
  53.  
  54. /* #define    MAX_ITER    10000 */
  55. static    int    max_iter = 10000;
  56. int    cg_num_iters;
  57.  
  58. /* matrix-as-routine type definition */
  59. /* #ifdef ANSI_C */
  60. /* typedef VEC    *(*MTX_FN)(void *params, VEC *x, VEC *out); */
  61. /* #else */
  62. typedef VEC    *(*MTX_FN)();
  63. /* #endif */
  64. #ifdef ANSI_C
  65. VEC    *spCHsolve(SPMAT *,VEC *,VEC *);
  66. #else
  67. VEC    *spCHsolve();
  68. #endif
  69.  
  70. /* cg_set_maxiter -- sets maximum number of iterations if numiter > 1
  71.     -- just returns current max_iter otherwise
  72.     -- returns old maximum */
  73. int    cg_set_maxiter(numiter)
  74. int    numiter;
  75. {
  76.     int    temp;
  77.  
  78.     if ( numiter < 2 )
  79.         return max_iter;
  80.     temp = max_iter;
  81.     max_iter = numiter;
  82.     return temp;
  83. }
  84.  
  85.  
  86. /* pccg -- solves A.x = b using pre-conditioner M
  87.             (assumed factored a la spCHfctr())
  88.     -- results are stored in x (if x != NULL), which is returned */
  89. VEC    *pccg(A,A_params,M_inv,M_params,b,eps,x)
  90. MTX_FN    A, M_inv;
  91. VEC    *b, *x;
  92. double    eps;
  93. void    *A_params, *M_params;
  94. {
  95.     VEC    *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  96.     int    k;
  97.     Real    alpha, beta, ip, old_ip, norm_b;
  98.  
  99.     if ( ! A || ! b )
  100.         error(E_NULL,"pccg");
  101.     if ( x == b )
  102.         error(E_INSITU,"pccg");
  103.     x = v_resize(x,b->dim);
  104.     if ( eps <= 0.0 )
  105.         eps = MACHEPS;
  106.  
  107.     r = v_get(b->dim);
  108.     p = v_get(b->dim);
  109.     q = v_get(b->dim);
  110.     z = v_get(b->dim);
  111.  
  112.     norm_b = v_norm2(b);
  113.  
  114.     v_zero(x);
  115.     r = v_copy(b,r);
  116.     old_ip = 0.0;
  117.     for ( k = 0; ; k++ )
  118.     {
  119.         if ( v_norm2(r) < eps*norm_b )
  120.             break;
  121.         if ( k > max_iter )
  122.             error(E_ITER,"pccg");
  123.         if ( M_inv )
  124.             (*M_inv)(M_params,r,z);
  125.         else
  126.             v_copy(r,z);    /* M == identity */
  127.         ip = in_prod(z,r);
  128.         if ( k )    /* if ( k > 0 ) ... */
  129.         {
  130.             beta = ip/old_ip;
  131.             p = v_mltadd(z,p,beta,p);
  132.         }
  133.         else        /* if ( k == 0 ) ... */
  134.         {
  135.             beta = 0.0;
  136.             p = v_copy(z,p);
  137.             old_ip = 0.0;
  138.         }
  139.         q = (*A)(A_params,p,q);
  140.         alpha = ip/in_prod(p,q);
  141.         x = v_mltadd(x,p,alpha,x);
  142.         r = v_mltadd(r,q,-alpha,r);
  143.         old_ip = ip;
  144.     }
  145.     cg_num_iters = k;
  146.  
  147.     V_FREE(p);
  148.     V_FREE(q);
  149.     V_FREE(r);
  150.     V_FREE(z);
  151.  
  152.     return x;
  153. }
  154.  
  155. /* sp_pccg -- a simple interface to pccg() which uses sparse matrix
  156.         data structures
  157.     -- assumes that LLT contains the Cholesky factorisation of the
  158.         actual pre-conditioner */
  159. VEC    *sp_pccg(A,LLT,b,eps,x)
  160. SPMAT    *A, *LLT;
  161. VEC    *b, *x;
  162. double    eps;
  163. {    return pccg(sp_mv_mlt,A,spCHsolve,LLT,b,eps,x);        }
  164.  
  165.  
  166. /*
  167.     Routines for performing the CGS (Conjugate Gradient Squared)
  168.     algorithm of P. Sonneveld:
  169.         "CGS, a fast Lanczos-type solver for nonsymmetric linear
  170.         systems", SIAM J. Sci. & Stat. Comp. v. 10, pp. 36--52
  171. */
  172.  
  173. /* cgs -- uses CGS to compute a solution x to A.x=b
  174.     -- the matrix A is not passed explicitly, rather a routine
  175.         A is passed where A(x,Ax,params) computes
  176.         Ax = A.x
  177.     -- the computed solution is passed */
  178. VEC    *cgs(A,A_params,b,r0,tol,x)
  179. MTX_FN    A;
  180. VEC    *x, *b;
  181. VEC    *r0;        /* tilde r0 parameter -- should be random??? */
  182. double    tol;        /* error tolerance used */
  183. void    *A_params;
  184. {
  185.     VEC    *p, *q, *r, *u, *v, *tmp1, *tmp2;
  186.     Real    alpha, beta, norm_b, rho, old_rho, sigma;
  187.     int    iter;
  188.  
  189.     if ( ! A || ! x || ! b || ! r0 )
  190.         error(E_NULL,"cgs");
  191.     if ( x->dim != b->dim || r0->dim != x->dim )
  192.         error(E_SIZES,"cgs");
  193.     if ( tol <= 0.0 )
  194.         tol = MACHEPS;
  195.  
  196.     p = v_get(x->dim);
  197.     q = v_get(x->dim);
  198.     r = v_get(x->dim);
  199.     u = v_get(x->dim);
  200.     v = v_get(x->dim);
  201.     tmp1 = v_get(x->dim);
  202.     tmp2 = v_get(x->dim);
  203.  
  204.     norm_b = v_norm2(b);
  205.     (*A)(A_params,x,tmp1);
  206.     v_sub(b,tmp1,r);
  207.     v_zero(p);    v_zero(q);
  208.     old_rho = 1.0;
  209.  
  210.     iter = 0;
  211.     while ( v_norm2(r) > tol*norm_b )
  212.     {
  213.         if ( ++iter > max_iter ) break;
  214.         /*    error(E_ITER,"cgs");  */
  215.         rho = in_prod(r0,r);
  216.         if ( old_rho == 0.0 )
  217.             error(E_SING,"cgs");
  218.         beta = rho/old_rho;
  219.         v_mltadd(r,q,beta,u);
  220.         v_mltadd(q,p,beta,tmp1);
  221.         v_mltadd(u,tmp1,beta,p);
  222.  
  223.         (*A)(A_params,p,v);
  224.  
  225.         sigma = in_prod(r0,v);
  226.         if ( sigma == 0.0 )
  227.             error(E_SING,"cgs");
  228.         alpha = rho/sigma;
  229.         v_mltadd(u,v,-alpha,q);
  230.         v_add(u,q,tmp1);
  231.  
  232.         (*A)(A_params,tmp1,tmp2);
  233.  
  234.         v_mltadd(r,tmp2,-alpha,r);
  235.         v_mltadd(x,tmp1,alpha,x);
  236.  
  237.         old_rho = rho;
  238.     }
  239.     cg_num_iters = iter;
  240.  
  241.     V_FREE(p);    V_FREE(q);    V_FREE(r);
  242.     V_FREE(u);    V_FREE(v);
  243.     V_FREE(tmp1);    V_FREE(tmp2);
  244.  
  245.     return x;
  246. }
  247.  
  248. /* sp_cgs -- simple interface for SPMAT data structures */
  249. VEC    *sp_cgs(A,b,r0,tol,x)
  250. SPMAT    *A;
  251. VEC    *b, *r0, *x;
  252. double    tol;
  253. {    return cgs(sp_mv_mlt,A,b,r0,tol,x);    }
  254.  
  255. /*
  256.     Routine for performing LSQR -- the least squares QR algorithm
  257.     of Paige and Saunders:
  258.         "LSQR: an algorithm for sparse linear equations and
  259.         sparse least squares", ACM Trans. Math. Soft., v. 8
  260.         pp. 43--71 (1982)
  261. */
  262. /* lsqr -- sparse CG-like least squares routine:
  263.     -- finds min_x ||A.x-b||_2 using A defined through A & AT
  264.     -- returns x (if x != NULL) */
  265. VEC    *lsqr(A,AT,A_params,b,tol,x)
  266. MTX_FN    A, AT;    /* AT is A transposed */
  267. VEC    *x, *b;
  268. double    tol;        /* error tolerance used */
  269. void    *A_params;
  270. {
  271.     VEC    *u, *v, *w, *tmp;
  272.     Real    alpha, beta, norm_b, phi, phi_bar,
  273.                 rho, rho_bar, rho_max, theta;
  274.     Real    s, c;    /* for Givens' rotations */
  275.     int    iter, m, n;
  276.  
  277.     if ( ! b || ! x )
  278.         error(E_NULL,"lsqr");
  279.     if ( tol <= 0.0 )
  280.         tol = MACHEPS;
  281.  
  282.     m = b->dim;    n = x->dim;
  283.     u = v_get((u_int)m);
  284.     v = v_get((u_int)n);
  285.     w = v_get((u_int)n);
  286.     tmp = v_get((u_int)n);
  287.     norm_b = v_norm2(b);
  288.  
  289.     v_zero(x);
  290.     beta = v_norm2(b);
  291.     if ( beta == 0.0 )
  292.         return x;
  293.     sv_mlt(1.0/beta,b,u);
  294.     tracecatch((*AT)(A_params,u,v),"lsqr");
  295.     alpha = v_norm2(v);
  296.     if ( alpha == 0.0 )
  297.         return x;
  298.     sv_mlt(1.0/alpha,v,v);
  299.     v_copy(v,w);
  300.     phi_bar = beta;        rho_bar = alpha;
  301.  
  302.     rho_max = 1.0;
  303.     iter = 0;
  304.     do {
  305.         if ( ++iter > max_iter )
  306.             error(E_ITER,"lsqr");
  307.  
  308.         tmp = v_resize(tmp,m);
  309.         tracecatch((*A) (A_params,v,tmp),"lsqr");
  310.  
  311.         v_mltadd(tmp,u,-alpha,u);
  312.         beta = v_norm2(u);    sv_mlt(1.0/beta,u,u);
  313.  
  314.         tmp = v_resize(tmp,n);
  315.         tracecatch((*AT)(A_params,u,tmp),"lsqr");
  316.         v_mltadd(tmp,v,-beta,v);
  317.         alpha = v_norm2(v);    sv_mlt(1.0/alpha,v,v);
  318.  
  319.         rho = sqrt(rho_bar*rho_bar+beta*beta);
  320.         if ( rho > rho_max )
  321.             rho_max = rho;
  322.         c   = rho_bar/rho;
  323.         s   = beta/rho;
  324.         theta   =  s*alpha;
  325.         rho_bar = -c*alpha;
  326.         phi     =  c*phi_bar;
  327.         phi_bar =  s*phi_bar;
  328.  
  329.         /* update x & w */
  330.         if ( rho == 0.0 )
  331.             error(E_SING,"lsqr");
  332.         v_mltadd(x,w,phi/rho,x);
  333.         v_mltadd(v,w,-theta/rho,w);
  334.     } while ( fabs(phi_bar*alpha*c) > tol*norm_b/rho_max );
  335.  
  336.     cg_num_iters = iter;
  337.  
  338.     V_FREE(tmp);    V_FREE(u);    V_FREE(v);    V_FREE(w);
  339.  
  340.     return x;
  341. }
  342.  
  343. /* sp_lsqr -- simple interface for SPMAT data structures */
  344. VEC    *sp_lsqr(A,b,tol,x)
  345. SPMAT    *A;
  346. VEC    *b, *x;
  347. double    tol;
  348. {    return lsqr(sp_mv_mlt,sp_vm_mlt,A,b,tol,x);    }
  349.  
  350.