home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / iter0.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  9KB  |  382 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.1 1994/01/13 05:38:23 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->A_par = NULL;
  107.    ip->ATx = ip->AT_par = NULL;
  108.    ip->Bx = ip->B_par = NULL;
  109.    ip->info = iter_std_info;
  110.    ip->stop_crit = iter_std_stop_crit;
  111.    ip->init_res = 0.0;
  112.    
  113.    return ip;
  114. }
  115.  
  116.  
  117. /* iter_free - release memory */
  118. int iter_free(ip)
  119. ITER *ip;
  120. {
  121.    if (ip == (ITER *)NULL) return -1;
  122.    
  123.    if (mem_info_is_on()) {
  124.       mem_bytes(TYPE_ITER,sizeof(ITER),0);
  125.       mem_numvar(TYPE_ITER,-1);
  126.    }
  127.  
  128.    if ( !ip->shared_x && ip->x != NULL ) v_free(ip->x);
  129.    if ( !ip->shared_b && ip->b != NULL ) v_free(ip->b);
  130.  
  131.    free((char *)ip);
  132.  
  133.    return 0;
  134. }
  135.  
  136. ITER *iter_resize(ip,new_lenb,new_lenx)
  137. ITER *ip;
  138. int new_lenb, new_lenx;
  139. {
  140.    VEC *old;
  141.  
  142.    if ( ip == (ITER *) NULL)
  143.      error(E_NULL,"iter_resize");
  144.  
  145.    old = ip->x;
  146.    ip->x = v_resize(ip->x,new_lenx);
  147.    if ( ip->shared_x && old != ip->x )
  148.      warning(WARN_SHARED_VEC,"iter_resize");
  149.    old = ip->b;
  150.    ip->b = v_resize(ip->b,new_lenb);
  151.    if ( ip->shared_b && old != ip->b )
  152.      warning(WARN_SHARED_VEC,"iter_resize");
  153.  
  154.    return ip;
  155. }
  156.  
  157.  
  158. /* print out ip structure - for diagnostic purposes mainly */
  159. void iter_dump(fp,ip)
  160. ITER *ip;
  161. FILE *fp;
  162. {
  163.    if (ip == NULL) {
  164.       fprintf(fp," ITER structure: NULL\n");
  165.       return;
  166.    }
  167.  
  168.    fprintf(fp,"\n ITER structure:\n");
  169.    fprintf(fp," ip->shared_x = %s, ip->shared_b = %s\n",
  170.        (ip->shared_x ? "TRUE" : "FALSE"),
  171.        (ip->shared_b ? "TRUE" : "FALSE") );
  172.    fprintf(fp," ip->k = %d, ip->limit = %d, ip->steps = %d, ip->eps = %g\n",
  173.        ip->k,ip->limit,ip->steps,ip->eps);
  174.    fprintf(fp," ip->x = 0x%p, ip->b = 0x%p\n",ip->x,ip->b);
  175.    fprintf(fp," ip->Ax = 0x%p, ip->A_par = 0x%p\n",ip->Ax,ip->A_par);
  176.    fprintf(fp," ip->ATx = 0x%p, ip->AT_par = 0x%p\n",ip->ATx,ip->AT_par);
  177.    fprintf(fp," ip->Bx = 0x%p, ip->B_par = 0x%p\n",ip->Bx,ip->B_par);
  178.    fprintf(fp," ip->info = 0x%p, ip->stop_crit = 0x%p, ip->init_res = %g\n",
  179.        ip->info,ip->stop_crit,ip->init_res);
  180.    fprintf(fp,"\n");
  181.    
  182. }
  183.  
  184.  
  185. /* copy the structure ip1 to ip2 preserving vectors x and b of ip2
  186.    (vectors x and b in ip2 are the same before and after iter_copy2)
  187.    if ip2 == NULL then a new structure is created with x and b being NULL
  188.    and other members are taken from ip1
  189. */
  190. ITER *iter_copy2(ip1,ip2)
  191. ITER *ip1, *ip2;
  192. {
  193.    VEC *x, *b;
  194.    int shx, shb;
  195.  
  196.    if (ip1 == (ITER *)NULL) 
  197.      error(E_NULL,"iter_copy2");
  198.  
  199.    if (ip2 == (ITER *)NULL) {
  200.       if ((ip2 = NEW(ITER)) == (ITER *) NULL)
  201.     error(E_MEM,"iter_copy2");
  202.       else if (mem_info_is_on()) {
  203.      mem_bytes(TYPE_ITER,0,sizeof(ITER));
  204.      mem_numvar(TYPE_ITER,1);
  205.       }
  206.       ip2->x = ip2->b = NULL;
  207.       ip2->shared_x = ip2->shared_x = FALSE;
  208.    }
  209.  
  210.    x = ip2->x;
  211.    b = ip2->b;
  212.    shb = ip2->shared_b;
  213.    shx = ip2->shared_x;
  214.    MEM_COPY(ip1,ip2,sizeof(ITER));
  215.    ip2->x = x;
  216.    ip2->b = b;
  217.    ip2->shared_x = shx;
  218.    ip2->shared_b = shb;
  219.  
  220.    return ip2;
  221. }
  222.  
  223.  
  224. /* copy the structure ip1 to ip2 copying also the vectors x and b */
  225. ITER *iter_copy(ip1,ip2)
  226. ITER *ip1, *ip2;
  227. {
  228.    VEC *x, *b;
  229.  
  230.    if (ip1 == (ITER *)NULL) 
  231.      error(E_NULL,"iter_copy");
  232.  
  233.    if (ip2 == (ITER *)NULL) {
  234.       if ((ip2 = NEW(ITER)) == (ITER *) NULL)
  235.     error(E_MEM,"iter_copy2");
  236.       else if (mem_info_is_on()) {
  237.      mem_bytes(TYPE_ITER,0,sizeof(ITER));
  238.      mem_numvar(TYPE_ITER,1);
  239.       }
  240.    }
  241.  
  242.    x = ip2->x;
  243.    b = ip2->b;
  244.  
  245.    MEM_COPY(ip1,ip2,sizeof(ITER));
  246.    if (ip1->x)
  247.      ip2->x = v_copy(ip1->x,x);
  248.    if (ip1->b)
  249.      ip2->b = v_copy(ip1->b,b);
  250.  
  251.    ip2->shared_x = ip2->shared_b = FALSE;
  252.  
  253.    return ip2;
  254. }
  255.  
  256.  
  257. /*** functions to generate sparse matrices with random entries ***/
  258.  
  259.  
  260. /* iter_gen_sym -- generate symmetric positive definite
  261.    n x n matrix, 
  262.    nrow - number of nonzero entries in a row
  263.    */
  264. SPMAT    *iter_gen_sym(n,nrow)
  265. int    n, nrow;
  266. {
  267.    SPMAT    *A;
  268.    VEC            *u;
  269.    Real       s1;
  270.    int        i, j, k, k_max;
  271.    
  272.    if (nrow <= 1) nrow = 2;
  273.    /* nrow should be even */
  274.    if ((nrow & 1)) nrow -= 1;
  275.    A = sp_get(n,n,nrow);
  276.    u = v_get(A->m);
  277.    v_zero(u);
  278.    for ( i = 0; i < A->m; i++ )
  279.    {
  280.       k_max = ((rand() >> 8) % (nrow/2));
  281.       for ( k = 0; k <= k_max; k++ )
  282.       {
  283.      j = (rand() >> 8) % A->n;
  284.      s1 = mrand();
  285.      sp_set_val(A,i,j,s1);
  286.      sp_set_val(A,j,i,s1);
  287.      u->ve[i] += fabs(s1);
  288.      u->ve[j] += fabs(s1);
  289.       }
  290.    }
  291.    /* ensure that A is positive definite */
  292.    for ( i = 0; i < A->m; i++ )
  293.      sp_set_val(A,i,i,u->ve[i] + 1.0);
  294.    
  295.    V_FREE(u);
  296.    return A;
  297. }
  298.  
  299.  
  300. /* iter_gen_nonsym -- generate non-symmetric m x n sparse matrix, m >= n 
  301.    nrow - number of entries in a row;
  302.    diag - number which is put in diagonal entries and then permuted
  303.    (if diag is zero then 1.0 is there)
  304. */
  305. SPMAT    *iter_gen_nonsym(m,n,nrow,diag)
  306. int    m, n, nrow;
  307. double diag;
  308. {
  309.    SPMAT    *A;
  310.    PERM        *px;
  311.    int        i, j, k, k_max;
  312.    Real        s1;
  313.    
  314.    if (nrow <= 1) nrow = 2;
  315.    if (diag == 0.0) diag = 1.0;
  316.    A = sp_get(m,n,nrow);
  317.    px = px_get(n);
  318.    for ( i = 0; i < A->m; i++ )
  319.    {
  320.       k_max = (rand() >> 8) % (nrow-1);
  321.       for ( k = 0; k <= k_max; k++ )
  322.       {
  323.      j = (rand() >> 8) % A->n;
  324.      s1 = mrand();
  325.      sp_set_val(A,i,j,-s1);
  326.       }
  327.    }
  328.    /* to make it likely that A is nonsingular, use pivot... */
  329.    for ( i = 0; i < 2*A->n; i++ )
  330.    {
  331.       j = (rand() >> 8) % A->n;
  332.       k = (rand() >> 8) % A->n;
  333.       px_transp(px,j,k);
  334.    }
  335.    for ( i = 0; i < A->n; i++ )
  336.      sp_set_val(A,i,px->pe[i],diag);  
  337.    
  338.    PX_FREE(px);
  339.    return A;
  340. }
  341.  
  342.  
  343. /* iter_gen_nonsym -- generate non-symmetric positive definite 
  344.    n x n sparse matrix;
  345.    nrow - number of entries in a row
  346. */
  347. SPMAT    *iter_gen_nonsym_posdef(n,nrow)
  348. int    n, nrow;
  349. {
  350.    SPMAT    *A;
  351.    PERM        *px;
  352.    VEC          *u;
  353.    int        i, j, k, k_max;
  354.    Real        s1;
  355.    
  356.    if (nrow <= 1) nrow = 2;
  357.    A = sp_get(n,n,nrow);
  358.    px = px_get(n);
  359.    u = v_get(A->m);
  360.    v_zero(u);
  361.    for ( i = 0; i < A->m; i++ )
  362.    {
  363.       k_max = (rand() >> 8) % (nrow-1);
  364.       for ( k = 0; k <= k_max; k++ )
  365.       {
  366.      j = (rand() >> 8) % A->n;
  367.      s1 = mrand();
  368.      sp_set_val(A,i,j,-s1);
  369.      u->ve[i] += fabs(s1);
  370.       }
  371.    }
  372.    /* ensure that A is positive definite */
  373.    for ( i = 0; i < A->m; i++ )
  374.      sp_set_val(A,i,i,u->ve[i] + 1.0);
  375.    
  376.    PX_FREE(px);
  377.    V_FREE(u);
  378.    return A;
  379. }
  380.  
  381.  
  382.