home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / iter0 < prev    next >
Encoding:
Text File  |  1994-08-03  |  8.8 KB  |  389 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. /* iter0.c  14/09/93 */
  28.  
  29. /* ITERATIVE METHODS - service functions */
  30.  
  31. /* functions for creating and releasing ITER structures;
  32.    for memory information;
  33.    for getting some values from an ITER variable;
  34.    for changing values in an ITER variable;
  35.    see also iter.c
  36. */
  37.  
  38. #include        <stdio.h>
  39. #include    <math.h>
  40. #include        "iter.h"
  41.  
  42.  
  43. static char rcsid[] = "$Id: iter0.c,v 1.2 1994/06/20 05:22:10 des Exp $";
  44.  
  45.  
  46. /* standard functions */
  47.  
  48. /* standard information */
  49. void iter_std_info(ip,nres,res,Bres)
  50. ITER *ip;
  51. double nres;
  52. VEC *res, *Bres;
  53. {
  54.    if (nres >= 0.0)
  55.      printf(" %d. residual = %g\n",ip->steps,nres);
  56.    else 
  57.      printf(" %d. residual = %g (WARNING !!! should be >= 0) \n",
  58.         ip->steps,nres);
  59. }
  60.  
  61. /* standard stopping criterion */
  62. int iter_std_stop_crit(ip, nres, res, Bres)
  63. ITER *ip;
  64. double nres;
  65. VEC *res, *Bres;
  66. {
  67.    /* save initial value of the residual in init_res */
  68.    if (ip->steps == 0)
  69.      ip->init_res = fabs(nres);
  70.  
  71.    /* standard stopping criterium */
  72.    if (nres <= ip->init_res*ip->eps) return TRUE; 
  73.    return FALSE;
  74. }
  75.  
  76.  
  77. /* iter_get - create a new structure pointing to ITER */
  78.  
  79. ITER *iter_get(lenb, lenx)
  80. int lenb, lenx;
  81. {
  82.    ITER *ip;
  83.  
  84.    if ((ip = NEW(ITER)) == (ITER *) NULL)
  85.      error(E_MEM,"iter_get");
  86.    else if (mem_info_is_on()) {
  87.       mem_bytes(TYPE_ITER,0,sizeof(ITER));
  88.       mem_numvar(TYPE_ITER,1);
  89.    }
  90.  
  91.    /* default values */
  92.    
  93.    ip->shared_x = FALSE;
  94.    ip->shared_b = FALSE;
  95.    ip->k = 0;
  96.    ip->limit = ITER_LIMIT_DEF;
  97.    ip->eps = ITER_EPS_DEF;
  98.    ip->steps = 0;
  99.  
  100.    if (lenb > 0) ip->b = v_get(lenb);
  101.    else ip->b = (VEC *)NULL;
  102.  
  103.    if (lenx > 0) ip->x = v_get(lenx);
  104.    else ip->x = (VEC *)NULL;
  105.  
  106.    ip->Ax = ip->ATx = ip->Bx = NULL;
  107.    ip->A_par = ip->AT_par = ip->B_par = NULL;
  108.    ip->info = iter_std_info;
  109.    ip->stop_crit = iter_std_stop_crit;
  110.    ip->init_res = 0.0;
  111.    
  112.    return ip;
  113. }
  114.  
  115.  
  116. /* iter_free - release memory */
  117. int iter_free(ip)
  118. ITER *ip;
  119. {
  120.    if (ip == (ITER *)NULL) return -1;
  121.    
  122.    if (mem_info_is_on()) {
  123.       mem_bytes(TYPE_ITER,sizeof(ITER),0);
  124.       mem_numvar(TYPE_ITER,-1);
  125.    }
  126.  
  127.    if ( !ip->shared_x && ip->x != NULL ) v_free(ip->x);
  128.    if ( !ip->shared_b && ip->b != NULL ) v_free(ip->b);
  129.  
  130.    free((char *)ip);
  131.  
  132.    return 0;
  133. }
  134.  
  135. ITER *iter_resize(ip,new_lenb,new_lenx)
  136. ITER *ip;
  137. int new_lenb, new_lenx;
  138. {
  139.    int old;
  140.  
  141.    if ( ip == (ITER *) NULL)
  142.      error(E_NULL,"iter_resize");
  143.  
  144.    if (new_lenx <= 0) ip->x = VNULL;
  145.    else {
  146.       old = ip->x->dim;
  147.       ip->x = v_resize(ip->x,new_lenx);
  148.       if ( ip->shared_x && old != new_lenx)
  149.     warning(WARN_SHARED_VEC,"iter_resize");
  150.    }
  151.    
  152.    if (new_lenb <= 0) ip->b = VNULL;
  153.    else {
  154.       old = ip->b->dim;
  155.       ip->b = v_resize(ip->b,new_lenb);
  156.       if ( ip->shared_b && old != new_lenb)
  157.     warning(WARN_SHARED_VEC,"iter_resize");
  158.    }
  159.       
  160.    return ip;
  161. }
  162.  
  163.  
  164.  
  165. /* print out ip structure - for diagnostic purposes mainly */
  166. void iter_dump(fp,ip)
  167. ITER *ip;
  168. FILE *fp;
  169. {
  170.    if (ip == NULL) {
  171.       fprintf(fp," ITER structure: NULL\n");
  172.       return;
  173.    }
  174.  
  175.    fprintf(fp,"\n ITER structure:\n");
  176.    fprintf(fp," ip->shared_x = %s, ip->shared_b = %s\n",
  177.        (ip->shared_x ? "TRUE" : "FALSE"),
  178.        (ip->shared_b ? "TRUE" : "FALSE") );
  179.    fprintf(fp," ip->k = %d, ip->limit = %d, ip->steps = %d, ip->eps = %g\n",
  180.        ip->k,ip->limit,ip->steps,ip->eps);
  181.    fprintf(fp," ip->x = 0x%p, ip->b = 0x%p\n",ip->x,ip->b);
  182.    fprintf(fp," ip->Ax = 0x%p, ip->A_par = 0x%p\n",ip->Ax,ip->A_par);
  183.    fprintf(fp," ip->ATx = 0x%p, ip->AT_par = 0x%p\n",ip->ATx,ip->AT_par);
  184.    fprintf(fp," ip->Bx = 0x%p, ip->B_par = 0x%p\n",ip->Bx,ip->B_par);
  185.    fprintf(fp," ip->info = 0x%p, ip->stop_crit = 0x%p, ip->init_res = %g\n",
  186.        ip->info,ip->stop_crit,ip->init_res);
  187.    fprintf(fp,"\n");
  188.    
  189. }
  190.  
  191.  
  192. /* copy the structure ip1 to ip2 preserving vectors x and b of ip2
  193.    (vectors x and b in ip2 are the same before and after iter_copy2)
  194.    if ip2 == NULL then a new structure is created with x and b being NULL
  195.    and other members are taken from ip1
  196. */
  197. ITER *iter_copy2(ip1,ip2)
  198. ITER *ip1, *ip2;
  199. {
  200.    VEC *x, *b;
  201.    int shx, shb;
  202.  
  203.    if (ip1 == (ITER *)NULL) 
  204.      error(E_NULL,"iter_copy2");
  205.  
  206.    if (ip2 == (ITER *)NULL) {
  207.       if ((ip2 = NEW(ITER)) == (ITER *) NULL)
  208.     error(E_MEM,"iter_copy2");
  209.       else if (mem_info_is_on()) {
  210.      mem_bytes(TYPE_ITER,0,sizeof(ITER));
  211.      mem_numvar(TYPE_ITER,1);
  212.       }
  213.       ip2->x = ip2->b = NULL;
  214.       ip2->shared_x = ip2->shared_x = FALSE;
  215.    }
  216.  
  217.    x = ip2->x;
  218.    b = ip2->b;
  219.    shb = ip2->shared_b;
  220.    shx = ip2->shared_x;
  221.    MEM_COPY(ip1,ip2,sizeof(ITER));
  222.    ip2->x = x;
  223.    ip2->b = b;
  224.    ip2->shared_x = shx;
  225.    ip2->shared_b = shb;
  226.  
  227.    return ip2;
  228. }
  229.  
  230.  
  231. /* copy the structure ip1 to ip2 copying also the vectors x and b */
  232. ITER *iter_copy(ip1,ip2)
  233. ITER *ip1, *ip2;
  234. {
  235.    VEC *x, *b;
  236.  
  237.    if (ip1 == (ITER *)NULL) 
  238.      error(E_NULL,"iter_copy");
  239.  
  240.    if (ip2 == (ITER *)NULL) {
  241.       if ((ip2 = NEW(ITER)) == (ITER *) NULL)
  242.     error(E_MEM,"iter_copy2");
  243.       else if (mem_info_is_on()) {
  244.      mem_bytes(TYPE_ITER,0,sizeof(ITER));
  245.      mem_numvar(TYPE_ITER,1);
  246.       }
  247.    }
  248.  
  249.    x = ip2->x;
  250.    b = ip2->b;
  251.  
  252.    MEM_COPY(ip1,ip2,sizeof(ITER));
  253.    if (ip1->x)
  254.      ip2->x = v_copy(ip1->x,x);
  255.    if (ip1->b)
  256.      ip2->b = v_copy(ip1->b,b);
  257.  
  258.    ip2->shared_x = ip2->shared_b = FALSE;
  259.  
  260.    return ip2;
  261. }
  262.  
  263.  
  264. /*** functions to generate sparse matrices with random entries ***/
  265.  
  266.  
  267. /* iter_gen_sym -- generate symmetric positive definite
  268.    n x n matrix, 
  269.    nrow - number of nonzero entries in a row
  270.    */
  271. SPMAT    *iter_gen_sym(n,nrow)
  272. int    n, nrow;
  273. {
  274.    SPMAT    *A;
  275.    VEC            *u;
  276.    Real       s1;
  277.    int        i, j, k, k_max;
  278.    
  279.    if (nrow <= 1) nrow = 2;
  280.    /* nrow should be even */
  281.    if ((nrow & 1)) nrow -= 1;
  282.    A = sp_get(n,n,nrow);
  283.    u = v_get(A->m);
  284.    v_zero(u);
  285.    for ( i = 0; i < A->m; i++ )
  286.    {
  287.       k_max = ((rand() >> 8) % (nrow/2));
  288.       for ( k = 0; k <= k_max; k++ )
  289.       {
  290.      j = (rand() >> 8) % A->n;
  291.      s1 = mrand();
  292.      sp_set_val(A,i,j,s1);
  293.      sp_set_val(A,j,i,s1);
  294.      u->ve[i] += fabs(s1);
  295.      u->ve[j] += fabs(s1);
  296.       }
  297.    }
  298.    /* ensure that A is positive definite */
  299.    for ( i = 0; i < A->m; i++ )
  300.      sp_set_val(A,i,i,u->ve[i] + 1.0);
  301.    
  302.    V_FREE(u);
  303.    return A;
  304. }
  305.  
  306.  
  307. /* iter_gen_nonsym -- generate non-symmetric m x n sparse matrix, m >= n 
  308.    nrow - number of entries in a row;
  309.    diag - number which is put in diagonal entries and then permuted
  310.    (if diag is zero then 1.0 is there)
  311. */
  312. SPMAT    *iter_gen_nonsym(m,n,nrow,diag)
  313. int    m, n, nrow;
  314. double diag;
  315. {
  316.    SPMAT    *A;
  317.    PERM        *px;
  318.    int        i, j, k, k_max;
  319.    Real        s1;
  320.    
  321.    if (nrow <= 1) nrow = 2;
  322.    if (diag == 0.0) diag = 1.0;
  323.    A = sp_get(m,n,nrow);
  324.    px = px_get(n);
  325.    for ( i = 0; i < A->m; i++ )
  326.    {
  327.       k_max = (rand() >> 8) % (nrow-1);
  328.       for ( k = 0; k <= k_max; k++ )
  329.       {
  330.      j = (rand() >> 8) % A->n;
  331.      s1 = mrand();
  332.      sp_set_val(A,i,j,-s1);
  333.       }
  334.    }
  335.    /* to make it likely that A is nonsingular, use pivot... */
  336.    for ( i = 0; i < 2*A->n; i++ )
  337.    {
  338.       j = (rand() >> 8) % A->n;
  339.       k = (rand() >> 8) % A->n;
  340.       px_transp(px,j,k);
  341.    }
  342.    for ( i = 0; i < A->n; i++ )
  343.      sp_set_val(A,i,px->pe[i],diag);  
  344.    
  345.    PX_FREE(px);
  346.    return A;
  347. }
  348.  
  349.  
  350. /* iter_gen_nonsym -- generate non-symmetric positive definite 
  351.    n x n sparse matrix;
  352.    nrow - number of entries in a row
  353. */
  354. SPMAT    *iter_gen_nonsym_posdef(n,nrow)
  355. int    n, nrow;
  356. {
  357.    SPMAT    *A;
  358.    PERM        *px;
  359.    VEC          *u;
  360.    int        i, j, k, k_max;
  361.    Real        s1;
  362.    
  363.    if (nrow <= 1) nrow = 2;
  364.    A = sp_get(n,n,nrow);
  365.    px = px_get(n);
  366.    u = v_get(A->m);
  367.    v_zero(u);
  368.    for ( i = 0; i < A->m; i++ )
  369.    {
  370.       k_max = (rand() >> 8) % (nrow-1);
  371.       for ( k = 0; k <= k_max; k++ )
  372.       {
  373.      j = (rand() >> 8) % A->n;
  374.      s1 = mrand();
  375.      sp_set_val(A,i,j,-s1);
  376.      u->ve[i] += fabs(s1);
  377.       }
  378.    }
  379.    /* ensure that A is positive definite */
  380.    for ( i = 0; i < A->m; i++ )
  381.      sp_set_val(A,i,i,u->ve[i] + 1.0);
  382.    
  383.    PX_FREE(px);
  384.    V_FREE(u);
  385.    return A;
  386. }
  387.  
  388.  
  389.