home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / lufactor < prev    next >
Encoding:
Text File  |  1994-03-13  |  6.4 KB  |  247 lines
  1. expt)
  2. VEC    *a;
  3. double    offset;
  4. int    *expt;
  5. {
  6.    Real    mant, tmp_fctr;
  7.    int    i, tmp_expt;
  8.    
  9.    if ( ! a )
  10.      error(E_NULL,"product");
  11.    
  12.    mant = 1.0;
  13.    *expt = 0;
  14.    if ( offset == 0.0 )
  15.      for ( i = 0; i < a->dim; i++ )
  16.      {
  17.     mant *= frexp(a->ve[i],&tmp_expt);
  18.     *expt += tmp_expt;
  19.     if ( ! (i % 10) )
  20.     {
  21.        mant = frexp(mant,&tmp_expt);
  22.        *expt += tmp_expt;
  23.     }
  24.      }
  25.    else
  26.      for ( i = 0; i < a->dim; i++ )
  27.      {
  28.     tmp_fctr = a->ve[i] - offset;
  29.     tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset :
  30.       MACHEPS*offset;
  31.     mant *= frexp(tmp_fctr,&tmp_expt);
  32.     *expt += tmp_expt;
  33.     if ( ! (i % 10) )
  34.     {
  35.        mant = frexp(mant,&tmp_expt);
  36.        *expt += tmp_expt;
  37.     }
  38.      }
  39.    
  40.    mant = frexp(mant,&tmp_expt);
  41.    *expt += tmp_expt;
  42.    
  43.    return mant;
  44. }
  45.  
  46. /* product2 -- returns the product of a long list of numbers
  47.    -- answer stored in mant (mantissa) and expt (exponent) */
  48. static    double    product2(a,k,expt)
  49. VEC    *a;
  50. int    k;    /* entry of a to leave out */
  51. int    *expt;
  52. {
  53.    Real    mant, mu, tmp_fctr;
  54.    int    i, tmp_expt;
  55.    
  56.    if ( ! a )
  57.      error(E_NULL,"product2");
  58.    if ( k < 0 || k >= a->dim )
  59.      error(E_BOUNDS,"product2");
  60.    
  61.    mant = 1.0;
  62.    *expt = 0;
  63.    mu = a->ve[k];
  64.    for ( i = 0; i < a->dim; i++ )
  65.    {
  66.       if ( i == k )
  67.     continue;
  68.       tmp_fctr = a->ve[i] - mu;
  69.       tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu;
  70.       mant *= frexp(tmp_fctr,&tmp_expt);
  71.       *expt += tmp_expt;
  72.       if ( ! (i % 10) )
  73.       {
  74.      mant = frexp(mant,&tmp_expt);
  75.      *expt += tmp_expt;
  76.       }
  77.    }
  78.    mant = frexp(mant,&tmp_expt);
  79.    *expt += tmp_expt;
  80.    
  81.    return mant;
  82. }
  83.  
  84. /* dbl_cmp -- comparison function to pass to qsort() */
  85. static    int    dbl_cmp(x,y)
  86. Real    *x, *y;
  87. {
  88.    Real    tmp;
  89.    
  90.    tmp = *x - *y;
  91.    return (tmp > 0 ? 1 : tmp < 0 ? -1: 0);
  92. }
  93.  
  94. /* iter_lanczos2 -- lanczos + error estimate for every e-val
  95.    -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978
  96.    -- returns multiple e-vals where multiple e-vals may not exist
  97.    -- returns evals vector */
  98. VEC    *iter_lanczos2(ip,evals,err_est)
  99. ITER     *ip;            /* ITER structure */
  100. VEC    *evals;        /* eigenvalue vector */
  101. VEC    *err_est;    /* error estimates of eigenvalues */
  102. {
  103.    VEC        *a;
  104.    static    VEC    *b=VNULL, *a2=VNULL, *b2=VNULL;
  105.    Real    beta, pb_mant, det_mant, det_mant1, det_mant2;
  106.    int    i, pb_expt, det_expt, det_expt1, det_expt2;
  107.    
  108.    if ( ! ip )
  109.      error(E_NULL,"iter_lanczos2");
  110.    if ( ! ip->Ax || ! ip->x )
  111.      error(E_NULL,"iter_lanczos2");
  112.    if ( ip->k <= 0 )
  113.      error(E_RANGE,"iter_lanczos2");
  114.    
  115.    a = evals;
  116.    a = v_resize(a,(u_int)ip->k);
  117.    b = v_resize(b,(u_int)(ip->k-1));
  118.    MEM_STAT_REG(b,TYPE_VEC);
  119.    
  120.    iter_lanczos(ip,a,b,&beta,MNULL);
  121.    
  122.    /* printf("# beta =%g\n",beta); */
  123.    pb_mant = 0.0;
  124.    if ( err_est )
  125.    {
  126.       pb_mant = product(b,(double)0.0,&pb_expt);
  127.       /* printf("# pb_mant = %g, pb_expt = %d\n",pb_mant, pb_expt); */
  128.    }
  129.    
  130.    /* printf("# diags =\n");    v_output(a); */
  131.    /* printf("# off diags =\n");    v_output(b); */
  132.    a2 = v_resize(a2,a->dim - 1);
  133.    b2 = v_resize(b2,b->dim - 1);
  134.    MEM_STAT_REG(a2,TYPE_VEC);
  135.    MEM_STAT_REG(b2,TYPE_VEC);
  136.    for ( i = 0; i < a2->dim - 1; i++ )
  137.    {
  138.       a2->ve[i] = a->ve[i+1];
  139.       b2->ve[i] = b->ve[i+1];
  140.    }
  141.    a2->ve[a2->dim-1] = a->ve[a2->dim];
  142.    
  143.    trieig(a,b,MNULL);
  144.    
  145.    /* sort evals as a courtesy */
  146.    qsort((void *)(a->ve),(int)(a->dim),sizeof(Real),(int (*)())dbl_cmp);
  147.    
  148.    /* error estimates */
  149.    if ( err_est )
  150.    {
  151.       err_est = v_resize(err_est,(u_int)ip->k);
  152.       
  153.       trieig(a2,b2,MNULL);
  154.       /* printf("# a =\n");    v_output(a); */
  155.       /* printf("# a2 =\n");    v_output(a2); */
  156.       
  157.       for ( i = 0; i < a->dim; i++ )
  158.       {
  159.      det_mant1 = product2(a,i,&det_expt1);
  160.      det_mant2 = product(a2,(double)a->ve[i],&det_expt2);
  161.      /* rights reserved.
  162. **
  163. **                 Meschach Library
  164. ** 
  165. ** This Meschach Library is provided "as is" without any express 
  166. ** or implied warranty of any kind with respect to this software. 
  167. ** In particular the authors shall not be liable for any direct, 
  168. ** indirect, special, incidental or consequential damages arising 
  169. ** in any way from use of the software.
  170. ** 
  171. ** Everyone is granted permission to copy, modify and redistribute this
  172. ** Meschach Library, provided:
  173. **  1.  All copies contain this copyright notice.
  174. **  2.  All modified copies shall carry a notice stating who
  175. **      made the last modification and the date of such modification.
  176. **  3.  No charge is made for this software or works derived from it.  
  177. **      This clause shall not be construed as constraining other software
  178. **      distributed on the same medium as this software, nor is a
  179. **      distribution fee considered a charge.
  180. **
  181. ***************************************************************************/
  182.  
  183.  
  184. /* matop.c 1.3 11/25/87 */
  185.  
  186.  
  187. #include    <stdio.h>
  188. #include    "matrix.h"
  189.  
  190. static    char    rcsid[] = "$Id: matop.c,v 1.3 1994/01/13 05:30:28 des Exp $";
  191.  
  192.  
  193. /* m_add -- matrix addition -- may be in-situ */
  194. MAT    *m_add(mat1,mat2,out)
  195. MAT    *mat1,*mat2,*out;
  196. {
  197.     u_int    m,n,i;
  198.  
  199.     if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  200.         error(E_NULL,"m_add");
  201.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  202.         error(E_SIZES,"m_add");
  203.     if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  204.         out = m_resize(out,mat1->m,mat1->n);
  205.     m = mat1->m;    n = mat1->n;
  206.     for ( i=0; i<m; i++ )
  207.     {
  208.         __add__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  209.         /**************************************************
  210.         for ( j=0; j<n; j++ )
  211.             out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
  212.         **************************************************/
  213.     }
  214.  
  215.     return (out);
  216. }
  217.  
  218. /* m_sub -- matrix subtraction -- may be in-situ */
  219. MAT    *m_sub(mat1,mat2,out)
  220. MAT    *mat1,*mat2,*out;
  221. {
  222.     u_int    m,n,i;
  223.  
  224.     if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  225.         error(E_NULL,"m_sub");
  226.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  227.     e for any direct, 
  228. ** indirect, special, incidental or consequential damages arising 
  229. ** in any way from use of the software.
  230. ** 
  231. ** Everyone is granted permission to copy, modify and redistribute this
  232. ** Meschach Library, provided:
  233. **  1.  All copies contain this copyright notice.
  234. **  2.  All modified copies shall carry a notice stating who
  235. **      made the last modification and the date of such modification.
  236. **  3.  No charge is made for this software or works derived from it.  
  237. **      This clause shall not be construed as constraining other software
  238. **      distributed on the same medium as this software, nor is a
  239. **      distribution fee considered a charge.
  240. **
  241. ***************************************************************************/
  242.  
  243.  
  244. /*
  245.         Files for matrix computations
  246.  
  247.