home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / itersym < prev    next >
Encoding:
Text File  |  1994-01-13  |  13.5 KB  |  619 lines
  1.  all rights reserved.
  2. **
  3. **                 Meschach Library
  4. ** 
  5. ** This Meschach Library is provided "as is" without any express 
  6. ** or implied warranty of any kind with respect to this software. 
  7. ** In particular the authors shall not be liable for any direct, 
  8. ** indirect, special, incidental or consequential damages arising 
  9. ** in any way from use of the software.
  10. ** 
  11. ** Everyone is granted permission to copy, modify and redistribute this
  12. ** Meschach Library, provided:
  13. **  1.  All copies contain this copyright notice.
  14. **  2.  All modified copies shall carry a notice stating who
  15. **      made the last modification and the date of such modification.
  16. **  3.  No charge is made for this software or works derived from it.  
  17. **      This clause shall not be construed as constraining other software
  18. **      distributed on the same medium as this software, nor is a
  19. **      distribution fee considered a charge.
  20. **
  21. ***************************************************************************/
  22.  
  23.  
  24. /*
  25.         Files for matrix computations
  26.  
  27.     Householder transformation file. Contains routines for calculating
  28.     householder transformations, applying them to vectors and matrices
  29.     by both row & column.
  30. */
  31.  
  32. /* hsehldr.c 1.3 10/8/87 */
  33. static    char    rcsid[] = "$Id: hsehldr.c,v 1.2 1994/01/13 05:36:29 des Exp $";
  34.  
  35. #include    <stdio.h>
  36. #include    <math.h>
  37. #include    "matrix.h"
  38. #include        "matrix2.h"
  39.  
  40.  
  41. /* hhvec -- calulates Householder vector to eliminate all entries after the
  42.     i0 entry of the vector vec. It is returned as out. May be in-situ */
  43. VEC    *hhvec(vec,i0,beta,out,newval)
  44. VEC    *vec,*out;
  45. u_int    i0;
  46. Real    *beta,*newval;
  47. {
  48.     Real    norm;
  49.  
  50.     out = _v_copy(vec,out,i0);
  51.     norm = sqrt(_in_prod(out,out,i0));
  52.     if ( norm <= 0.0 )
  53.     {
  54.         *beta = 0.0;
  55.         return (out);
  56.     }
  57.     *beta = 1.0/(norm * (norm+fabs(out->ve[i0])));
  58.     if ( out->ve[i0] > 0.0 )
  59.         *newval = -norm;
  60.     else
  61.         *newval = norm;
  62.     out->ve[i0] -= *newval;
  63.  
  64.     return (out);
  65. }
  66.  
  67. /* hhtrvec -- apply Householder transformation to vector -- may be in-situ */
  68. VEC    *hhtrvec(hh,beta,i0,in,out)
  69. VEC    *hh,*in,*out;    /* hh = Householder vector */
  70. u_int    i0;
  71. double    beta;
  72. {
  73.     Real    scale;
  74.     /* u_int    i; */
  75.  
  76.     if ( hh==(VEC *)NULL || in==(VEC *)NULL )
  77.         error(E_NULL,"hhtrvec");
  78.     if ( in->dim != hh->dim )
  79.         error(E_SIZES,"hhtrvec");
  80.     if ( i0 > in->dim )
  81.         error(E_BOUNDS,"hhtrvec");
  82.  
  83.     scale = beta*_in_prod(hh,in,i0);
  84.     out = v_copy(in,out);
  85.     __mltadd__(&(out->ve[i0]),&(hh->ve[i0]),-scale,(int)(in->dim-i0));
  86.     /************************************************************
  87.     for ( i=i0; i<in->dim; i++ )
  88.         out->ve[i] = in->ve[i] - scale*hh->ve[i];
  89.     ************************************************************/
  90.  
  91.     return (out);
  92. }
  93.  
  94. /* hhtrrows -- transform a matrix by a Householder vector by rows
  95.     starting at row i0 from column j0 -- in-situ */
  96. MAT    *hhtrrows(M,i0,j0,hh,beta)
  97. MAT    *M;
  98. u_int    i0, j0;
  99. VEC    *hh;
  100. double    beta;
  101. {
  102.     Real    ip, scale;
  103.     int    i /*, j */;
  104.  
  105.     if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  106.         error(E_NULL,"hhtrrows");
  107.     if ( M->n != hh->dim )
  108.         error(E_RANGE,"hhtrrows");
  109.     if ( i0 > M->m || j0 > M->n )
  110.         error(E_BOUNDS,"hhtrrows");
  111.  
  112.     if ( beta == 0.0 )    return (M);
  113.  
  114.     /* for each row ... */
  115.     for ( i = i0; i < M->m; i++ )
  116.     {    /* compute inner product */
  117.         ip = __ip__(&(M->me[i][j0]),&(hh->ve[j0]),(int)(M->n-j0));
  118.         /**************************************************
  119.         ip = 0.0;
  120.         for ( j = j0; j < M->n; j++ )
  121.             ip += M->me[i][j]*hh->ve[j];
  122.         **************************************************/
  123.         scale = beta*ip;
  124.         if ( scale == 0.0 )
  125.             continue;
  126.  
  127.         /* do operation */
  128.         __mltadd__(&(M->me[i][j0]),&(hh->ve[j0]),-scale,
  129.                             (int)(M->n-j0));
  130.         /**************************************************
  131.         for ( j = j0; j < M->n; j++ )
  132.             M->me[i][j] -= scale*hh->ve[j];
  133.         **************************************************/
  134.     }
  135.  
  136.     return (M);
  137. }
  138.  
  139.  
  140. /* hhtrcols -- transform a matrix by a Householder vector by columns
  141.     starting at row i0 from column j0 -- in-situ */
  142. MAT    *hhtrcols(M,i0,j0,hh,beta)
  143. MAT    *M;
  144. u_int    i0, j0;
  145. VEC    *hh;
  146. double    beta;
  147. {
  148.     /* Real    ip, scale; */
  149.     int    i /*, k */;
  150.     static    VEC    *w = VNULL;
  151.  
  152.     if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  153.         error(E_NULL,"hhtrcols");
  154.     if ( M->m != hh->dim )
  155.         error(E_SIZES,"hhtrcols");
  156.     if ( i0 > M->m || j0 > M->n )
  157.         error(E_BOUNDS,"hhtrcols");
  158.  
  159.     if ( beta == 0.0 )    return (M);
  160.  
  161.     w = v_resize(w,M->n);
  162.     MEM_STAT_REG(w,TYPE_VEC);
  163.     v_zero(w);
  164.  
  165.     for ( i = i0; i < M->m; i++ )
  166.         if ( hh->ve[i] != 0.0 )
  167.         __mltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
  168.                             (int)(M->n-j0));
  169.     for ( i = i0; i < M->m; i++ )
  170.         if ( hh->ve[i] != 0.0 )
  171.         __mltadd__(&(M->me[i][j0]),&(w->ve[j0]),-beta*hh->ve[i],
  172.                             (int)(M->n-j0));
  173.     return (M);
  174. }
  175.  
  176. FileDataŵinitEÿÿÿè%@Ðë
  177. /**************************************************************************
  178. **
  179. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  180. **
  181. **                 Meschach Library
  182. ** 
  183. ** This Meschach Library is provided "as is" without any express 
  184. ** or implied warranty of any kind with respect to this software. 
  185. ** In particular the authors shall not be liable for any direct, 
  186. ** indirect, special, incidental or consequential damages arising 
  187. ** in any way from use of the software.
  188. ** 
  189. ** Everyone is granted permission to copy, modify and redistribute this
  190. ** Meschach Library, provided:
  191. **  1.  All copies contain this copyright notice.
  192. **  2.  All modified copies shall carry a notice stating who
  193. **      made the last modification and the date of such modification.
  194. **  3.  No charge is made for this software or works derived from it.  
  195. **      This clause shall not be construed as constraining other software
  196. **      distributed on the same medium as this software, nor is a
  197. **      distribution fee considered a charge.
  198. **
  199. ***************************************************************************/
  200.  
  201.  
  202. /*
  203.     This is a file of routines for zero-ing, and initialising
  204.     vectors, matrices and permutations.
  205.     This is to be included in the matrix.a library
  206. */
  207.  
  208. static    char    rcsid[] = "$Id: init.c,v 1.6 1994/01/13 05:36:58 des Exp $";
  209.  
  210. #include    <stdio.h>
  211. #include    "matrix.h"
  212.  
  213. /* v_zero -- zero the vector x */
  214. VEC    *v_zero(x)
  215. VEC    *x;
  216. {
  217.     if ( x == VNULL )
  218.         error(E_NULL,"v_zero");
  219.  
  220.     __zero__(x->ve,x->dim);
  221.     /* for ( i = 0; i < x->dim; i++ )
  222.         x->ve[i] = 0.0; */
  223.  
  224.     return x;
  225. }
  226.  
  227.  
  228. /* iv_zero -- zero the vector ix */
  229. IVEC    *iv_zero(ix)
  230. IVEC    *ix;
  231. {
  232.    int i;
  233.    
  234.    if ( ix == IVNULL )
  235.      error(E_NULL,"iv_zero");
  236.    
  237.    for ( i = 0; i < ix->dim; i++ )
  238.      ix->ive[i] = 0; 
  239.    
  240.    return ix;
  241. }
  242.  
  243.  
  244. /* m_zero -- zero the matrix A */
  245. MAT    *m_zero(A)
  246. MAT    *A;
  247. {
  248.     int    i, A_m, A_n;
  249.     Real    **A_me;
  250.  
  251.     if ( A == MNULL )
  252.         error(E_NULL,"m_zero");
  253.  
  254.     A_m = A->m;    A_n = A->n;    A_me = A->me;
  255.     for ( i = 0; i < A_m; i++ )
  256.         __zero__(A_me[i],A_n);
  257.         /* for ( j = 0; j < A_n; j++ )
  258.             A_me[i][j] = 0.0; */
  259.  
  260.     return A;
  261. }
  262.  
  263. /* mat_id -- set A to being closest to identity matrix as possible
  264.     -- i.e. A[i][j] == 1 if i == j and 0 otherwise */
  265. MAT    *m_ident(A)
  266. MAT    *A;
  267. {
  268.     int    i, size;
  269.  
  270.     if ( A == MNULL )
  271.         error(E_NULL,"m_ident");
  272.  
  273.     m_zero(A);
  274.     size = min(A->m,A->n);
  275.     for ( i = 0; i < size; i++ )
  276.         A->me[i][i] = 1.0;
  277.  
  278.     return A;
  279. }
  280.     
  281. /* px_ident -- set px to identity permutation */
  282. PERM    *px_ident(px)
  283. PERM    *px;
  284. {
  285.     int    i, px_size;
  286.     u_int    *px_pe;
  287.  
  288.     if ( px == PNULL )
  289.         error(E_NULL,"px_ident");
  290.  
  291.     px_size = px->size;    px_pe = px->pe;
  292.     for ( i = 0; i < px_size; i++ )
  293.         px_pe[i] = i;
  294.  
  295.     return px;
  296. }
  297.  
  298. /* Pseudo random number generator data structures */
  299. /* Knuth's lagged Fibonacci-based generator: See "Seminumerical Algorithms:
  300.    The Art of Computer Programming" sections 3.2-3.3 */
  301.  
  302. #ifdef ANSI_C
  303. #ifndef LONG_MAX
  304. #include    <limits.h>
  305. #endif
  306. #endif
  307.  
  308. #ifdef LONG_MAX
  309. #define MODULUS    LONG_MAX
  310. #else
  311. #define MODULUS    1000000000L    /* assuming long's at least 32 bits long */
  312. #endif
  313. #define MZ    0L
  314.  
  315. static long mrand_list[56];
  316. static int  started = FALSE;
  317. static int  inext = 0, inextp = 31;
  318.  
  319.  
  320. /* mrand -- pseudo-random number generator */
  321. #ifdef ANSI_C
  322. double mrand(void)
  323. #else
  324. double mrand()
  325. #endif
  326. {
  327.     long    lval;
  328.     static Real  factor = 1.0/((Real)MODULUS);
  329.     
  330.     if ( ! started )
  331.     smrand(3127);
  332.     
  333.     inext = (inext >= 54) ? 0 : inext+1;
  334.     inextp = (inextp >= 54) ? 0 : inextp+1;
  335.  
  336.     lval = mrand_list[inext]-mrand_list[inextp];
  337.     if ( lval < 0L )
  338.     lval += MODULUS;
  339.     mrand_list[inext] = lval;
  340.     
  341.     return (double)lval*factor;
  342. }
  343.  
  344. /* mrandlist -- fills the array a[] with len random numbers */
  345. void    mrandlist(a, len)
  346. Real    a[];
  347. int    len;
  348. {
  349.     int        i;
  350.     long    lval;
  351.     static Real  factor = 1.0/((Real)MODULUS);
  352.     
  353.     if ( ! started )
  354.     smrand(3127);
  355.     
  356.     for ( i = 0; i < len; i++ )
  357.     {
  358.     inext = (inext >= 54) ? 0 : inext+1;
  359.     inextp = (inextp >= 54) ? 0 : inextp+1;
  360.     
  361.     lval = mrand_list[inext]-mrand_list[inextp];
  362.     if ( lval < 0L )
  363.         lval += MODULUS;
  364.     mrand_list[inext] = lval;
  365.     
  366.     a[i] = (Real)lval*factor;
  367.     }
  368. }
  369.  
  370.  
  371. /* smrand -- set seed for mrand() */
  372. void smrand(seed)
  373. int    seed;
  374. {
  375.     int        i;
  376.  
  377.     mrand_list[0] = (123413*seed) % MODULUS;
  378.     for ( i = 1; i < 55; i++ )
  379.     mrand_list[i] = (123413*mrand_list[i-1]) % MODULUS;
  380.  
  381.     started = TRUE;
  382.  
  383.     /* run mrand() through the list sufficient times to
  384.        thoroughly randomise the array */
  385.     for ( i = 0; i < 55*55; i++ )
  386.     mrand();
  387. }
  388. #undef MODULUS
  389. #undef MZ
  390. #undef FAC
  391.  
  392. /* v_rand -- initialises x to be a random vector, components
  393.     independently & uniformly ditributed between 0 and 1 */
  394. VEC    *v_rand(x)
  395. VEC    *x;
  396. {
  397.     /* int    i; */
  398.  
  399.     if ( ! x )
  400.         error(E_NULL,"v_rand");
  401.  
  402.     /* for ( i = 0; i < x->dim; i++ ) */
  403.         /* x->ve[i] = rand()/((Real)MAX_RAND); */
  404.         /* x->ve[i] = mrand(); */
  405.     mrandlist(x->ve,x->dim);
  406.  
  407.     return x;
  408. }
  409.  
  410. /* m_rand -- initialises A to be a random vector, components
  411.     independently & uniformly distributed between 0 and 1 */
  412. MAT    *m_rand(A)
  413. MAT    *A;
  414. {
  415.     int    i /* , j */;
  416.  
  417.     if ( ! A )
  418.         error(E_NULL,"m_rand");
  419.  
  420.     for ( i = 0; i < A->m; i++ )
  421.         /* for ( j = 0; j < A->n; j++ ) */
  422.             /* A->me[i][j] = rand()/((Real)MAX_RAND); */
  423.             /* A->me[i][j] = mrand(); */
  424.         mrandlist(A->me[i],A->n);
  425.  
  426.     return A;
  427. }
  428.  
  429. /* v_ones -- fills x with one's */
  430. VEC    *v_ones(x)
  431. VEC    *x;
  432. {
  433.     int    i;
  434.  
  435.     if ( ! x )
  436.         error(E_NULL,"v_ones");
  437.  
  438.     for ( i = 0; i < x->dim; i++ )
  439.         x->ve[i] = 1.0;
  440.  
  441.     return x;
  442. }
  443.  
  444. /* m_ones -- fills matrix with one's */
  445. MAT    *m_ones(A)
  446. MAT    *A;
  447. {
  448.     int    i, j;
  449.  
  450.     if ( ! A )
  451.         error(E_NULL,"m_ones");
  452.  
  453.     for ( i = 0; i < A->m; i++ )
  454.         for ( j = 0; j < A->n; j++ )
  455.             A->me[i][j] = 1.0;
  456.  
  457.     return A;
  458. }
  459.  
  460. /* v_count -- initialises x so that x->ve[i] == i */
  461. VEC    *v_count(x)
  462. VEC    *x;
  463. {
  464.     int    i;
  465.  
  466.     if ( ! x )
  467.         error(E_NULL,"v_count");
  468.  
  469.     for ( i = 0; i < x->d printf("# det_mant1=%g, det_expt1=%d\n",
  470.         det_mant1,det_expt1); */
  471.      /* printf("# det_mant2=%g, det_expt2=%d\n",
  472.         det_mant2,det_expt2); */
  473.      if ( det_mant1 == 0.0 )
  474.      {   /* multiple e-val of T */
  475.         err_est->ve[i] = 0.0;
  476.         continue;
  477.      }
  478.      else if ( det_mant2 == 0.0 )
  479.      {
  480.         err_est->ve[i] = HUGE;
  481.         continue;
  482.      }
  483.      if ( (det_expt1 + det_expt2) % 2 )
  484.        /* if odd... */
  485.        det_mant = sqrt(2.0*fabs(det_mant1*det_mant2));
  486.      else /* if even... */
  487.        det_mant = sqrt(fabs(det_mant1*det_mant2));
  488.      det_expt = (det_expt1+det_expt2)/2;
  489.      err_est->ve[i] = fabs(beta*
  490.                    ldexp(pb_mant/det_mant,pb_expt-det_expt));
  491.       }
  492.    }
  493.    
  494.    return a;
  495. }
  496.  
  497. /* iter_splanczos2 -- version of iter_lanczos2() that uses sparse matrix data
  498.    structure */
  499.  
  500. VEC    *iter_splanczos2(A,m,x0,evals,err_est)
  501. SPMAT    *A;
  502. int     m;
  503. VEC    *x0;        /* initial vector */
  504. VEC    *evals;        /* eigenvalue vector */
  505. VEC    *err_est;    /* error estimates of eigenvalues */
  506. {    
  507.    ITER *ip;
  508.    VEC *a;
  509.    
  510.    ip = iter_get(0,0);
  511.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  512.    ip->A_par = (void *) A;
  513.    ip->x = x0;
  514.    ip->k = m;
  515.    a = iter_lanczos2(ip,evals,err_est);    
  516.    ip->shared_x = ip->shared_b = TRUE;
  517.    iter_free(ip);   /* release only ITER structure */
  518.    return a;
  519. }
  520.  
  521.  
  522.  
  523.  
  524. /*
  525.   Conjugate gradient method
  526.   Another variant - mainly for testing
  527.   */
  528.  
  529. VEC  *iter_cg1(ip)
  530. ITER *ip;
  531. {
  532.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  533.    Real    alpha;
  534.    double inner,nres;
  535.    VEC *rr;   /* rr == r or rr == z */
  536.    
  537.    if (ip == INULL)
  538.      error(E_NULL,"iter_cg");
  539.    if (!ip->Ax || !ip->b)
  540.      error(E_NULL,"iter_cg");
  541.    if ( ip->x == ip->b )
  542.      error(E_INSITU,"iter_cg");
  543.    if (!ip->stop_crit)
  544.      error(E_NULL,"iter_cg");
  545.    
  546.    if ( ip->eps <= 0.0 )
  547.      ip->eps = MACHEPS;
  548.    
  549.    r = v_resize(r,ip->b->dim);
  550.    p = v_resize(p,ip->b->dim);
  551.    q = v_resize(q,ip->b->dim);
  552.    
  553.    MEM_STAT_REG(r,TYPE_VEC);
  554.    MEM_STAT_REG(p,TYPE_VEC);
  555.    MEM_STAT_REG(q,TYPE_VEC);
  556.    
  557.    if (ip->Bx != (Fun_Ax)NULL) {
  558.       z = v_resize(z,ip->b->dim);
  559.       MEM_STAT_REG(z,TYPE_VEC);
  560.       rr = z;
  561.    }
  562.    else rr = r;
  563.    
  564.    if (ip->x != VNULL) {
  565.       if (ip->x->dim != ip->b->dim)
  566.     error(E_SIZES,"iter_cg");
  567.       ip->Ax(ip->A_par,ip->x,p);            /* p = A*x */
  568.       v_sub(ip->b,p,r);                 /* r = b - A*x */
  569.    }
  570.    else {  /* ip->x == 0 */
  571.       ip->x = v_get(ip->b->dim);
  572.       ip->shared_x = FALSE;
  573.       v_copy(ip->b,r);
  574.    }
  575.    
  576.    if (ip->Bx) (ip->Bx)(ip->B_par,r,p);
  577.    else v_copy(r,p);
  578.    
  579.    inner = in_prod(p,r);
  580.    nres = sqrt(fabs(inner));
  581.    if (ip->info) ip->info(ip,nres,r,p);
  582.    if ( ip->stop_crit(ip,nres,r,p) ) return ip->x;
  583.    
  584.    for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
  585.    {
  586.       ip->Ax(ip->A_par,p,q);
  587.       inner = in_prod(q,p);
  588.       if (inner <= 0.0) {
  589.      warning(WARN_RES_LESS_0,"iter_cg");
  590.      break;
  591.       }
  592.       alpha = in_prod(p,r)/inner;
  593.       v_mltadd(ip->x,p,alpha,ip->x);
  594.       v_mltadd(r,q,-alpha,r);
  595.       
  596.       rr = r;
  597.       if (ip->Bx) {
  598.      ip->Bx(ip->B_par,r,z);
  599.      rr = z;
  600.       }
  601.       
  602.       nres = in_prod(r,rr);
  603.       if (nres < 0.0) {
  604.      warning(WARN_RES_LESS_0,"iter_cg");
  605.      break;
  606.       }
  607.       nres = sqrt(fabs(nres));
  608.       if (ip->info) ip->info(ip,nres,r,z);
  609.       if ( ip->stop_crit(ip,nres,r,z) ) break;
  610.       
  611.       alpha = -in_prod(rr,q)/inner;
  612.       v_mltadd(rr,p,alpha,p);
  613.       
  614.    }
  615.    
  616.    return ip->x;
  617. }
  618.  
  619.