home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / ivecop < prev    next >
Encoding:
Text File  |  1994-01-13  |  9.6 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. /* ivecop.c  */
  28.  
  29. #include    <stdio.h>
  30. #include     "matrix.h"
  31.  
  32. static    char    rcsid[] = "$Id: ivecop.c,v 1.5 1994/01/13 05:45:30 des Exp $";
  33.  
  34. static char    line[MAXLINE];
  35.  
  36.  
  37.  
  38. /* iv_get -- get integer vector -- see also memory.c */
  39. IVEC    *iv_get(dim)
  40. int    dim;
  41. {
  42.    IVEC    *iv;
  43.    /* u_int    i; */
  44.    
  45.    if (dim < 0)
  46.      error(E_NEG,"iv_get");
  47.  
  48.    if ((iv=NEW(IVEC)) == IVNULL )
  49.      error(E_MEM,"iv_get");
  50.    else if (mem_info_is_on()) {
  51.       mem_bytes(TYPE_IVEC,0,sizeof(IVEC));
  52.       mem_numvar(TYPE_IVEC,1);
  53.    }
  54.    
  55.    iv->dim = iv->max_dim = dim;
  56.    if ((iv->ive = NEW_A(dim,int)) == (int *)NULL )
  57.      error(E_MEM,"iv_get");
  58.    else if (mem_info_is_on()) {
  59.       mem_bytes(TYPE_IVEC,0,dim*sizeof(int));
  60.    }
  61.    
  62.    return (iv);
  63. }
  64.  
  65. /* iv_free -- returns iv & asoociated memory back to memory heap */
  66. int    iv_free(iv)
  67. IVEC    *iv;
  68. {
  69.    if ( iv==IVNULL || iv->dim > MAXDIM )
  70.      /* don't trust it */
  71.      return (-1);
  72.    
  73.    if ( iv->ive == (int *)NULL ) {
  74.       if (mem_info_is_on()) {
  75.      mem_bytes(TYPE_IVEC,sizeof(IVEC),0);
  76.      mem_numvar(TYPE_IVEC,-1);
  77.       }
  78.       free((char *)iv);
  79.    }
  80.    else
  81.    {
  82.       if (mem_info_is_on()) {
  83.      mem_bytes(TYPE_IVEC,sizeof(IVEC)+iv->max_dim*sizeof(int),0);
  84.      mem_numvar(TYPE_IVEC,-1);
  85.       }    
  86.       free((char *)iv->ive);
  87.       free((char *)iv);
  88.    }
  89.    
  90.    return (0);
  91. }
  92.  
  93. /* iv_resize -- returns the IVEC with dimension new_dim
  94.    -- iv is set to the zero vector */
  95. IVEC    *iv_resize(iv,new_dim)
  96. IVEC    *iv;
  97. int    new_dim;
  98. {
  99.    int    i;
  100.    
  101.    if (new_dim < 0)
  102.      error(E_NEG,"iv_resize");
  103.  
  104.    if ( ! iv )
  105.      return iv_get(new_dim);
  106.    
  107.    if (new_dim == iv->dim)
  108.      return iv;
  109.  
  110.    if ( new_dim > iv->max_dim )
  111.    {
  112.       if (mem_info_is_on()) {
  113.      mem_bytes(TYPE_IVEC,iv->max_dim*sizeof(int),
  114.               new_dim*sizeof(int));
  115.       }
  116.       iv->ive = RENEW(iv->ive,new_dim,int);
  117.       if ( ! iv->ive )
  118.     error(E_MEM,"iv_resize");
  119.       iv->max_dim = new_dim;
  120.    }
  121.    if ( iv->dim <= new_dim )
  122.      for ( i = iv->dim; i < new_dim; i++ )
  123.        iv->ive[i] = 0;
  124.    iv->dim = new_dim;
  125.    
  126.    return iv;
  127. }
  128.  
  129. /* iv_copy -- copy integer vector in to out
  130.    -- out created/resized if necessary */
  131. IVEC    *iv_copy(in,out)
  132. IVEC    *in, *out;
  133. {
  134.    int        i;
  135.    
  136.    if ( ! in )
  137.      error(E_NULL,"iv_copy");
  138.    out = iv_resize(out,in->dim);
  139.    for ( i = 0; i < in->dim; i++ )
  140.      out->ive[i] = in->ive[i];
  141.    
  142.  Also the member ip->init_res is updated indirectly by 
  143.      ip->stop_crit.
  144.      */
  145.       if (ip->steps == 0) {                /* information for a user */
  146.      if (ip->info) (*ip->info)(ip,nres,As,rr); 
  147.      if ( (*ip->stop_crit)(ip,nres,As,rr) ) { 
  148.           /* iterative process is finished */
  149.         done = TRUE; 
  150.         break;
  151.      }
  152.       }
  153.       else if (nres <= 0.0) 
  154.     break;  /* residual is zero -> finish the process */ 
  155.       
  156.       /* save this residual in the first row of N */
  157.       v.ve = N->me[0];
  158.       v_copy(rr,&v);
  159.       
  160.       for (i = 0; i < ip->k && ip->steps <= ip->limit; i++) {
  161.      ip->steps++;
  162.      v.ve = N->me[i];                /* pointer to a row of N (=s_i) */
  163.      /* note that we must use here &v, not v */
  164.      (*ip->Ax)(ip->A_par,&v,As); 
  165.      rr = As;                        /* As = A*s_i */
  166.      if (ip->Bx) {
  167.         (*ip->Bx)(ip->B_par,As,z);    /* z = B*A*s_i  */
  168.         rr = z;
  169.      }
  170.      
  171.      if (i < ip->k - 1) {
  172.         s.ve = N->me[i+1];         /* pointer to a row of N (=s_{i+1}) */
  173.         v_copy(rr,&s);                   /* s_{i+1} = B*A*s_i */
  174.         for (j = 0; j <= i-1; j++) {
  175.            v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
  176.            beta->ve[j] = in_prod(&v,rr);    /* beta_{j,i} */
  177.                    /* s_{i+1} -= beta_{j,i}*s_{j+1} */
  178.            v_mltadd(&s,&v,- beta->ve[j],&s);    
  179.         }
  180.         
  181.          /* beta_{i,i} = ||s_{i+1}||_2 */
  182.         beta->ve[i] = nres = v_norm2(&s);     
  183.         if ( nres <= 0.0) break;         /* if s_{i+1} == 0 then finish */
  184.         sv_mlt(1.0/nres,&s,&s);           /* normalize s_{i+1} */
  185.         
  186.         v.ve = N->me[0];
  187.         alpha->ve[i] = in_prod(&v,&s);     /* alpha_i = (s_0 , s_{i+1}) */
  188.         
  189.      }
  190.      else {
  191.         for (j = 0; j <= i-1; j++) {
  192.            v.ve = N->me[j+1];      /* pointer to a row of N (=s_{j+1}) */
  193.            beta->ve[j] = in_prod(&v,rr);       /* beta_{j,i} */
  194.         }
  195.         
  196.         nres = in_prod(rr,rr);                 /* rr = B*A*s_{k-1} */
  197.         for (j = 0; j <= i-1; j++)
  198.               nres -= beta->ve[j]*beta->ve[j];
  199.         if (nres <= 0.0)  break;               /* s_k is zero */
  200.         else beta->ve[i] = sqrt(nres);         /* beta_{k-1,k-1} */
  201.         
  202.         v.ve = N->me[0];
  203.         alpha->ve[i] = in_prod(&v,rr); 
  204.         for (j = 0; j <= i-1; j++)
  205.               alpha->ve[i] -= beta->ve[j]*alpha->ve[j];
  206.         alpha->ve[i] /= beta->ve[i];                /* alpha_{k-1} */
  207.         
  208.      }
  209.      
  210.      set_col(H,i,beta);
  211.      
  212.      dd -= alpha->ve[i]*alpha->ve[i];
  213.      nres = sqrt(fabs((double) dd));
  214.      if (dd < 0.0)  {     /* do restart */
  215.         if (ip->info) (*ip->info)(ip,-nres,VNULL,VNULL);  
  216.         break;
  217.      }
  218.      
  219.      if (ip->info) (*ip->info)(ip,nres,VNULL,VNULL);     
  220.      if ( ip->stop_crit(ip,nres,VNULL,VNULL) ) {
  221.         /* stopping criterion is satisfied */
  222.         done = TRUE;
  223.         break;
  224.      }
  225.      
  226.       } /* end of for */
  227.       
  228.       if (nres <= 0.0) {
  229.      i--;
  230.      done = TRUE;
  231.       }
  232.       if (i >= ip->k) i = ip->k - 1;
  233.       
  234.       /* use (i+1) by (i+1) submatrix of H */
  235.       H = m_resize(H,i+1,i+1);
  236.       alpha = v_resize(alpha,i+1);
  237.       Usolve(H,alpha,alpha,0.0);       /* c_i is saved in alpha */
  238.       
  239.       for (j = 0; j <= i; j++) {
  240.      v.ve = N->me[j];
  241.      v_mltadd(ip->x,&v,alpha->ve[j],ip->x);
  242.       }
  243.       
  244.       
  245.       if (done) break;              /* stop the iterative process */
  246.       alpha = v_resize(alpha,ip->k);
  247.       H = m_resize(H,ip->k,ip->k);
  248.       
  249.    }  /* end of while */
  250.    
  251.    return ip->x;                    /* return the solution */
  252. }
  253.  
  254.  
  255.  
  256. /* iter_spmgcr - a simple interface to iter_mgcr */
  257. /* no preconditioner */
  258. VEC    *iter_spmgcr(A,B,b,tol,x,k,limit,steps)
  259. SPMAT    *A, *B;
  260. VEC    *b, *x;
  261. double    tol;
  262. int *steps,k,limit;
  263. {
  264.    ITER *ip;
  265.    
  266.    ip = iter_get(0,0);
  267.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  268.    ip->A_par = (void *) A;
  269.    if (B) {
  270.       ip->Bx = (Fun_Ax) sp_mv_mlt;
  271.       ip->B_par = (void *) B;
  272.    }
  273.    else {
  274.       ip->Bx = (Fun_Ax) NULL;
  275.       ip->B_par = NULL;
  276.    }
  277.  
  278.    ip->k = k;
  279.    ip->limit = limit;
  280.    ip->info = (Fun_info) NULL;
  281.    ip->b = b;
  282.    ip->eps = tol;
  283.    ip->x = x;
  284.    iter_mgcr(ip);
  285.    x = ip->x;
  286.    if (steps) *steps = ip->steps;
  287.    ip->shared_x = ip->shared_b = TRUE;
  288.    iter_free(ip);   /* release only ITER structure */
  289.    return x;        
  290. }
  291.  
  292.  
  293.  
  294. /* 
  295.   Conjugate gradients method for a normal equation
  296.   a preconditioner B must be symmetric !!
  297. */
  298. VEC  *iter_cgne(ip)
  299. ITER *ip;
  300. {
  301.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  302.    Real    alpha, beta, inner, old_inner, nres;
  303.    VEC *rr1;   /* pointer only */
  304.    
  305.    if (ip == INULL)
  306.      error(E_NULL,"iter_cgne");
  307.    if (!ip->Ax || ! ip->ATx || !ip->b)
  308.      error(E_NULL,"iter_cgne");
  309.    if ( ip->x == ip->b )
  310.      error(E_INSITU,"iter_cgne");
  311.    if (!ip->stop_crit)
  312.      error(E_NULL,"iter_cgne");
  313.    
  314.    if ( ip->eps <= 0.0 )
  315.      ip->eps = MACHEPS;
  316.    
  317.    r = v_resize(r,ip->b->dim);
  318.    p = v_resize(p,ip->b->dim);
  319.    q = v_resize(q,ip->b->dim);
  320.  
  321.    MEM_STAT_REG(r,TYPE_VEC);
  322.    MEM_STAT_REG(p,TYPE_VEC);
  323.    MEM_STAT_REG(q,TYPE_VEC);
  324.  
  325.    z = v_resize(z,ip->b->dim);
  326.    MEM_STAT_REG(z,TYPE_VEC);
  327.  
  328.    if (ip->x) {
  329.       if (ip->x->dim != ip->b->dim)
  330.     error(E_SIZES,"iter_cgne");
  331.       ip->Ax(ip->A_par,ip->x,p);            /* p = A*x */
  332.       v_sub(ip->b,p,z);                 /* z = b - A*x */
  333.    }
  334.    else {  /* ip->x == 0 */
  335.       ip->x = v_get(ip->b->dim);
  336.       ip->shared_x = FALSE;
  337.       v_copy(ip->b,z);
  338.    }
  339.    rr1 = z;
  340.    if (ip->Bx) {
  341.       (ip->Bx)(ip->B_par,rr1,p);
  342.       rr1 = p;
  343.    }
  344.    (ip->ATx)(ip->AT_par,rr1,r);        /* r = A^T*B*(b-A*x)  */
  345.  
  346.  
  347.    old_inner = 0.0;
  348.    for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
  349.    {
  350.       rr1 = r;
  351.       if ( ip->Bx ) {
  352.      (ip->Bx)(ip->B_par,r,z);        /* rr = B*r */
  353.      rr1 = z;
  354.       }
  355.  
  356.       inner = in_prod(r,rr1);
  357.       nres = sqrt(fabs(inner));
  358.       if (ip->info) ip->info(ip,nres,r,rr1);
  359.       if ( ip->stop_crit(ip,nres,r,rr1) ) break;
  360.  
  361.       if ( ip->steps )    /* if ( ip->steps > 0 ) ... */
  362.       {
  363.      beta = inner/old_inner;
  364.      p = v_mltadd(rr1,p,beta,p);
  365.       }
  366.       else        /* if ( ip->steps == 0 ) ... */
  367.       {
  368.      beta = 0.0;
  369.      p = v_copy(rr1,p);
  370.      old_inner = 0.0;
  371.       }
  372.       (ip->Ax)(ip->A_par,p,q);     /* q = A*p */
  373.       if (ip->Bx) {
  374.      (ip->Bx)(ip->B_par,q,z);
  375.      (ip->ATx)(ip->AT_par,z,q);
  376.      rr1 = q;            /* q = A^T*B*A*p */
  377.       }
  378.       else {
  379.      (ip->ATx)(ip->AT_par,q,z);    /* z = A^T*A*p */
  380.      rr1 = z;
  381.       }
  382.  
  383.       alpha = inner/in_prod(rr1,p);
  384.       v_mltadd(ip->x,p,alpha,ip->x);
  385.       v_mltadd(r,rr1,-alpha,r);
  386.       old_inner = inner;
  387.    }
  388.  
  389.    return ip->