home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / lanczos < prev    next >
Encoding:
Text File  |  1994-01-13  |  7.6 KB  |  324 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.     File containing Lanczos type routines for finding eigenvalues
  29.     of large, sparse, symmetic matrices
  30. */
  31.  
  32. #include    <stdio.h>
  33. #include    <math.h>
  34. #include    "matrix.h"
  35. #include    "sparse.h"
  36.  
  37. static char rcsid[] = "$Id: lanczos.c,v 1.4 1994/01/13 05:28:24 des Exp $";
  38.  
  39. #ifdef ANSI_C
  40. extern    VEC    *trieig(VEC *,VEC *,MAT *);
  41. #else
  42. extern    VEC    *trieig();
  43. #endif
  44.  
  45. /* lanczos -- raw lanczos algorithm -- no re-orthogonalisation
  46.     -- creates T matrix of size == m,
  47.         but no larger than before beta_k == 0
  48.     -- uses passed routine to do matrix-vector multiplies */
  49. void    lanczos(A_fn,A_params,m,x0,a,b,beta2,Q)
  50. VEC    *(*A_fn)();    /* VEC *(*A_fn)(void *A_params,VEC *in, VEC *out) */
  51. void    *A_params;
  52. int    m;
  53. VEC    *x0, *a, *b;
  54. Real    *beta2;
  55. MAT    *Q;
  56. {
  57.     int    j;
  58.     VEC    *v, *w, *tmp;
  59.     Real    alpha, beta;
  60.  
  61.     if ( ! A_fn || ! x0 || ! a || ! b )
  62.         error(E_NULL,"lanczos");
  63.     if ( m <= 0 )
  64.         error(E_BOUNDS,"lanczos");
  65.     if ( Q && ( Q->m < x0->dim || Q->n < m ) )
  66.         error(E_SIZES,"lanczos");
  67.  
  68.     a = v_resize(a,(u_int)m);    b = v_resize(b,(u_int)(m-1));
  69.     v = v_get(x0->dim);
  70.     w = v_get(x0->dim);
  71.     tmp = v_get(x0->dim);
  72.  
  73.     beta = 1.0;
  74.     /* normalise x0 as w */
  75.     sv_mlt(1.0/v_norm2(x0),x0,w);
  76.  
  77.     (*A_fn)(A_params,w,v);
  78.  
  79.     for ( j = 0; j < m; j++ )
  80.     {
  81.         /* store w in Q if Q not NULL */
  82.         if ( Q )
  83.             set_col(Q,j,w);
  84.  
  85.         alpha = in_prod(w,v);
  86.         a->ve[j] = alpha;
  87.         v_mltadd(v,w,-alpha,v);
  88.         beta = v_norm2(v);
  89.         if ( beta == 0.0 )
  90.         {
  91.             v_resize(a,(u_int)j+1);
  92.             v_resize(b,(u_int)j);
  93.             *beta2 = 0.0;
  94.             if ( Q )
  95.             Q = m_resize(Q,Q->m,j+1);
  96.             return;
  97.         }
  98.         if ( j < m-1 )
  99.             b->ve[j] = beta;
  100.         v_copy(w,tmp);
  101.         sv_mlt(1/beta,v,w);
  102.         sv_mlt(-beta,tmp,v);
  103.         (*A_fn)(A_params,w,tmp);
  104.         v_add(v,tmp,v);
  105.     }
  106.     *beta2 = beta;
  107.  
  108.  
  109.     V_FREE(v);    V_FREE(w);    V_FREE(tmp);
  110. }
  111.  
  112. extern    double    frexp(), ldexp();
  113.  
  114. /* product -- returns the product of a long list of numbers
  115.     -- answer stored in mant (mantissa) and expt (exponent) */
  116. static    double    product(a,offset,expt)
  117. VEC    *a;
  118. double    offset;
  119. int    *expt;
  120. {
  121.     Real    mant, tmp_fctr;
  122.     int    i, tmp_expt;
  123.  
  124.     if ( ! a )
  125.         error(E_NULL,"product");
  126.  
  127.     mant = 1.0;
  128.     *expt = 0;
  129.     if ( offset == 0.0 )
  130.         for ( i = 0; i < a->dim; i++ )
  131.         {
  132.             mant *= frexp(a->ve[i],&tmp_expt);
  133.             *expt += tmp_expt;
  134.             if ( ! (i % 10) )
  135.             {
  136.                 mant = frexp(mant,&tmp_expt);
  137.                 *expt += tmp_expt;
  138.             }
  139.         }
  140.     else
  141.         for ( i = 0; i < a->dim; i++ )
  142.         {
  143.             tmp_fctr = a->ve[i] - offset;
  144.             tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset :
  145.                              MACHEPS*offset;
  146.             mant *= frexp(tmp_fctr,&tmp_expt);
  147.             *expt += tmp_expt;
  148.             if ( ! (i % 10) )
  149.             {
  150.                 mant = frexp(mant,&tmp_expt);
  151.                 *expt += tmp_expt;
  152.             }
  153.         }
  154.  
  155.     mant = frexp(mant,&tmp_expt);
  156.     *expt += tmp_expt;
  157.  
  158.     return mant;
  159. }
  160.  
  161. /* product2 -- returns the product of a long list of numbers
  162.     -- answer stored in mant (mantissa) and expt (exponent) */
  163. static    double    product2(a,k,expt)
  164. VEC    *a;
  165. int    k;    /* entry of a to leave out */
  166. int    *expt;
  167. {
  168.     Real    mant, mu, tmp_fctr;
  169.     int    i, tmp_expt;
  170.  
  171.     if ( ! a )
  172.         error(E_NULL,"product2");
  173.     if ( k < 0 || k >= a->dim )
  174.         error(E_BOUNDS,"product2");
  175.  
  176.     mant = 1.0;
  177.     *expt = 0;
  178.     mu = a->ve[k];
  179.     for ( i = 0; i < a->dim; i++ )
  180.     {
  181.         if ( i == k )
  182.             continue;
  183.         tmp_fctr = a->ve[i] - mu;
  184.         tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu;
  185.         mant *= frexp(tmp_fctr,&tmp_expt);
  186.         *expt += tmp_expt;
  187.         if ( ! (i % 10) )
  188.         {
  189.             mant = frexp(mant,&tmp_expt);
  190.             *expt += tmp_expt;
  191.         }
  192.     }
  193.     mant = frexp(mant,&tmp_expt);
  194.     *expt += tmp_expt;
  195.  
  196.     return mant;
  197. }
  198.  
  199. /* dbl_cmp -- comparison function to pass to qsort() */
  200. static    int    dbl_cmp(x,y)
  201. Real    *x, *y;
  202. {
  203.     Real    tmp;
  204.  
  205.     tmp = *x - *y;
  206.     return (tmp > 0 ? 1 : tmp < 0 ? -1: 0);
  207. }
  208.  
  209. /* lanczos2 -- lanczos + error estimate for every e-val
  210.     -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978
  211.     -- returns multiple e-vals where multiple e-vals may not exist
  212.     -- returns evals vector */
  213. VEC    *lanczos2(A_fn,A_params,m,x0,evals,err_est)
  214. VEC    *(*A_fn)();
  215. void    *A_params;
  216. int    m;
  217. VEC    *x0;        /* initial vector */
  218. VEC    *evals;        /* eigenvalue vector */
  219. VEC    *err_est;    /* error estimates of eigenvalues */
  220. {
  221.     VEC        *a;
  222.     static    VEC    *b=VNULL, *a2=VNULL, *b2=VNULL;
  223.     Real    beta, pb_mant, det_mant, det_mant1, det_mant2;
  224.     int    i, pb_expt, det_expt, det_expt1, det_expt2;
  225.  
  226.     if ( ! A_fn || ! x0 )
  227.         error(E_NULL,"lanczos2");
  228.     if ( m <= 0 )
  229.         error(E_RANGE,"lanczos2");
  230.  
  231.     a = evals;
  232.     a = v_resize(a,(u_int)m);
  233.     b = v_resize(b,(u_int)(m-1));
  234.     MEM_STAT_REG(b,TYPE_VEC);
  235.  
  236.     lanczos(A_fn,A_params,m,x0,a,b,&beta,MNULL);
  237.  
  238.     /* printf("# beta =%g\n",beta); */
  239.     pb_mant = 0.0;
  240.     if ( err_est )
  241.     {
  242.         pb_mant = product(b,(double)0.0,&pb_expt);
  243.         /* printf("# pb_mant = %g, pb_expt = %d\n",pb_mant, pb_expt); */
  244.     }
  245.  
  246.     /* printf("# diags =\n");    out_vec(a); */
  247.     /* printf("# off diags =\n");    out_vec(b); */
  248.     a2 = v_resize(a2,a->dim - 1);
  249.     b2 = v_resize(b2,b->dim - 1);
  250.     MEM_STAT_REG(a2,TYPE_VEC);
  251.     MEM_STAT_REG(b2,TYPE_VEC);
  252.     for ( i = 0; i < a2->dim - 1; i++ )
  253.     {
  254.         a2->ve[i] = a->ve[i+1];
  255.         b2->ve[i] = b->ve[i+1];
  256.     }
  257.     a2->ve[a2->dim-1] = a->ve[a2->dim];
  258.  
  259.     trieig(a,b,MNULL);
  260.  
  261.     /* sort evals as a courtesy */
  262.     qsort((void *)(a->ve),(int)(a->dim),sizeof(Real),(int (*)())dbl_cmp);
  263.  
  264.     /* error estimates */
  265.     if ( err_est )
  266.     {
  267.         err_est = v_resize(err_est,(u_int)m);
  268.  
  269.         trieig(a2,b2,MNULL);
  270.         /* printf("# a =\n");    out_vec(a); */
  271.         /* printf("# a2 =\n");    out_vec(a2); */
  272.  
  273.         for ( i = 0; i < a->dim; i++ )
  274.         {
  275.             det_mant1 = product2(a,i,&det_expt1);
  276.             det_mant2 = product(a2,(double)a->ve[i],&det_expt2);
  277.             /* printf("# det_mant1=%g, det_expt1=%d\n",
  278.                     det_mant1,det_expt1); */
  279.             /* printf("# det_mant2=%g, det_expt2=%d\n",
  280.                     det_mant2,det_expt2); */
  281.             if ( det_mant1 == 0.0 )
  282.             {   /* multiple e-val of T */
  283.                 err_est->ve[i] = 0.0;
  284.                 continue;
  285.             }
  286.             else if ( det_mant2 == 0.0 )
  287.             {
  288.                 err_est->ve[i] = HUGE;
  289.                 continue;
  290.             }
  291.             if ( (det_expt1 + det_expt2) % 2 )
  292.                 /* if odd... */
  293.                 det_mant = sqrt(2.0*fabs(det_mant1*det_mant2));
  294.             else /* if even... */
  295.                 det_mant = sqrt(fabs(det_mant1*det_mant2));
  296.             det_expt = (det_expt1+det_expt2)/2;
  297.             err_est->ve[i] = fabs(beta*
  298.                 ldexp(pb_mant/det_mant,pb_expt-det_expt));
  299.         }
  300.     }
  301.  
  302.     return a;
  303. }
  304.  
  305. /* sp_lanczos -- version that uses sparse matrix data structure */
  306. void    sp_lanczos(A,m,x0,a,b,beta2,Q)
  307. SPMAT    *A;
  308. int     m;
  309. VEC     *x0, *a, *b;
  310. Real  *beta2;
  311. MAT     *Q;
  312. {    lanczos(sp_mv_mlt,A,m,x0,a,b,beta2,Q);    }
  313.  
  314. /* sp_lanczos2 -- version of lanczos2() that uses sparse matrix data
  315.                     structure */
  316. VEC    *sp_lanczos2(A,m,x0,evals,err_est)
  317. SPMAT    *A;
  318. int    m;
  319. VEC    *x0;        /* initial vector */
  320. VEC    *evals;        /* eigenvalue vector */
  321. VEC    *err_est;    /* error estimates of eigenvalues */
  322. {    return lanczos2(sp_mv_mlt,A,m,x0,evals,err_est);    }
  323.  
  324.